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COMPUTER  SOLUTIONS  FOR  THERMAL-ACOUSTICAL 
OSCILLATIONS  IN  GAS-FILLED  TUBES 

M.    T.    Norton  and  R.    C.    Muhlenhaupt 


A  digital  computer  program  to  determine  solutions  for  thermal- 
acoustical  oscillations  in  gas -filled  pipes  is  described.      Using  a  typical 
temperature  gradient  and  tube  length  from  test  data,    the  program  cal- 
culates the  effect  of  changes  in  heat  transfer  coefficient,    friction  factor, 
and  pipe  diameter.     Details  of  the  program  are  explained.     A  comparison 
is  made  between  computation  and  one  particular  test  datum  point.     Al- 
though the  calculations  are  based  on  helium  as  the  media,    the  program 
can  accommodate  any  fluid  treated  as  a  perfect  gas. 

Key  Words:     Cryogenics,    thermal  oscillations,    liquid  helium 

1.     Introduction 
When  the  open  end  of  a  pipe  is  at  low  temperature  and  the  closed 
end  is  much  warmer,    thermal-acoustical  oscillations  often  occur.      These 
oscillations  are  common  in  cryogenic  vent  lines,    pressure  sensing  lines, 
and  liquid  transfer  lines.     This  study  is  an  investigation  of  this  usually 
undesirable  phenomenon.     The  result  is  a  digital  computer  program  to 
predict  frequencies  and  overpressure  amplitudes  from  pipe  diameter, 
temperature  gradient,    and  gas  properties.     One  test  datum  point  is  used 
to  establish  the  relation  between  text  book  friction  factors  and  friction 
level  in  an  oscillating  pipe. 

Curves  are  presented  showing  the  effect  of  pipe  diameter,    friction 
factor,    and  heat  transfer  coefficient  on  the  magnitude  of  the  oscillations. 

Although  the  program  is  specifically  for  helium  gas,    minor  changes 
can  adapt  it  to  other  gases. 

The  calculations  for  frequency  and  amplitude  are  performed  by 
quite  different  methods.     The  frequency  depends  on  acoustic  velocity, 


and  acoustic  frequencies  in  general  are  essentially  independent  of  ampli- 
tude.    So  it  is  assumed  that  in  thermal  oscillations,    the  frequency  is  inde- 
pendent of  amplitude  and,for  a  given  gas, depends  only  on  the  temperature 
gradient.     For  convenience,    when  determining  frequency,    an  initial  pulse 
of  pressure  is  assumed. 

The  amplitude  depends  on  the  frequency,    but  the  method  of  calcu- 
lating is  not  the  same  as  that  for  frequency.     In  amplitude  calculations 
various  assumptions  are  made  in  order  to  approximate  the  shape  of  the 
pressure  profile. 

2.     General  Theory  of  Amplitude  Calculations 

There  are  infinite  varieties  of  duct  configurations  which  can  pro- 
duce oscillations  when  inserted  into  a  cryogenic  dewar.     Of  these  the 
simplest  and  most  common  case  is  a  constant  diameter  tube,    closed  at 
the  warm  end,    and  open  at  the  cold  end.      This  is  the  same  as  a  dewar 
vent  line  with  the  shut  off  valve  closed,    or  a  pressure  measurement  line. 
It  is  well  known  that  these  are  subject  to  thermal  oscillations.      The 
general  theory  applicable  to  this  case  is  applicable  to  all  cases. 

The  problem  is  not  new  nor  is  it  confined  to  cryogenics.     Glass- 
blowers  frequently  experience  it.     A  description  of  the  phenomenon  and 
a  qualitative  explanation  is  found  on  page  230  of  Rayleigh  (1945). 

A  qualitative  explanation  of  thermal  oscillations,    fundamentally  in 
agreement  with  Rayleigh,    follows.     During  the  pressure  rise  portion  of 
an  oscillation  cycle  there  is  movement  of  gas  from  the  open  towards  the 
closed  end.     As  a  gas  element  travels,    its  pressure  increases  and,    from 
adiabatic  compression,    its  temperature  rises.      The  wall  temperature 
increases  toward  the  closed  end,    and  the  gradient  is  steeper  than  the 
temperature  gradient  in  the  gas.      Therefore,    there  is  a  temperature 
difference  tending  to  transfer  heat  from  the  wall  to  the  gas.     Because  heat 


transfer  requires  time,    first  the  pressure  increases  and  then  heat  is 
transferred.     After  the  pressure  in  the  element  has  peaked,    the  process 
is  reversed.      Then  there  is  a  series  of  pressure  drops  and  heat  transfers 
out  of  the  gas.     When  plotted  on  temperature-entropy  coordinates  this 
sequence  of  events  has  positive  work  area,    a  simplified  version  of  which 
is  shown  on  figure  1.      Thus  a  gas  element  in  contact  with  a  wall,   whose 
temperature  is  greater  at  the  closed  end,    and  whose  gradient  is   steeper 
than  the  adiabatic  gas  temperature  rise,    does  work  on  its  surroundings. 
If  each  gas  element  is  in  contact  with  this  type  of  temperature  gradient 
then,    clearly,   work  is  done  by  the  gas  elements  in  the  pipe.     Conversely, 
any  gas  element  in  contact  with  a  wall,    whose  temperature  gradient  is 
reversed  or  is  less  steep  than  adiabatic  compression,   will  have  work  done 
on  it.     Furthermore,    if  all  gas  elements  are  in  contact  with  a  gradient 
such  that  work  is  done  by  all  elements,    it  is  clear  that  the  motive  power 
exists  for  continued  oscillations.      Likewise,    if  all  gas  elements  are  in 
contact  with  a  gradient  such  that  work  is  absorbed  by  all  elements,    it  is 
clear  that  any  accidental  oscillations  will  quickly  damp  out.     For  the 
intermediate  case,   where  the  wall  temperature  gradient  sometimes  causes 
work  to  be  done  by,  and  sometimes  on, the  various  gas  elements,    it  is  not 
clear  if  chance  oscillations  damp  out  or  continue.     A  digital  computer 
program  can  be  used  to  sum  up  all  the  infinitesimals  of  work  done  in  one 
cycle  in  all  the  gas  elements  in  the  pipe. 

Thus,    in  an  oscillation  cycle,    positive  work  is  a  necessary  criterion 
for  the  existence  ol  thermal  oscillations.     If  the  work  summation  is  posi- 
tive,   any  accidental  oscillations  build  up  to  an  equilibrium  steady  state 
condition  or  until  heat  transfer  from  the  duct  wall  alters  the  temperature 
gradient,    in  which  case  the  oscillations  may  cease.     If  the  work  summation 
is  negative,    any  accidental  oscillations  quickly  damp  out. 


kj 


ENTROPY 


Figure  1.      Work  Cycle  Of  Simplified  Thermal  Oscillation  Cycle 


In  case  the  work  summation  is  positive,    at  equilibrium  it  must  ex- 
actly balance  losses.     There  are  two  types  of  losses,    friction  and  kinetic 
energy,    expelled  into  the  dewar.      Part  of  the  friction  loss  is  simply  wall 
friction  opposing  the  gas  motion  in  either  direction.     Also,   when  the  gas 
enters  the  pipe  from  the  dewar  there  is  turbulence  created  which  is  equiv- 
alent to  additional  wall  friction.     In  this  study  the  entrance  loss  will  be 
handled  as  an  additive  to  the  wall  friction. 

In  one  pressure  cycle  each  element  will  undergo  heat  transfer  and 
will  do  some  positive  or  negative  work.     As  each  element  slides  along  the 
wall  or  into  the  duct,    there  is  frictional  return  of  work  to  heat  in  the  gas 
element.     After  summing  up  for  all  elements  in  one  cycle,    the  excess  of 
positive  work  done  over  frictional  work  must  equal  the  kinetic  energy 
ejected  into  the  dewar  during  one  cycle. 

This  criterion,    that  the  net  work  done  after  subtracting  friction 
must  equal  the  kinetic  energy  ejected,    is  used  in  an  iterative  calculation 
to  determine  the  pressure  amplitude  of  helium  gas  in  a  given  pipe  with  a 
given  temperature  gradient.      The  method  used  is  to  estimate  the  pressure 
amplitude  and  calculate  for  all  gas  elements  in  the  duct,  the  work  done 
and  the  friction  and  kinetic  energy  ejected.     If  the  net  work  exceeds  the 
kinetic  energy,    the  estimate  is  low.  If  the  net  work  is  less  than  the  kinetic 
energy,    the  estimate  is  high.     With  a  few  iterations  it  is  possible  to  deter- 
mine the  pressure  amplitude  for  which  the  net  work  equals  the  kinetic 
energy  ejected.      This  is  the  calculated  pressure. 

Figure  2  illustrates  the  method  used.    The  particular  pipe  diameter 
under  investigation  is  .3175  cm  (1/8  inch).     In  this  illustration,    pipe  fric- 
tion is  calculated  using  friction  factors     that  remain  constant  during  each 
run,    and  heat  transfer  is  computed  using  unmodified  heat  transfer  coeffi- 
cients from  McAdams   (1954).     When  the  overpressure  estimate  is  low, 
the  ratio  of  ejected  kinetic  energy  to  work  done  is  less  than  unity.     When 


»- 

o 

ii    a  t 
»~  o   E 

.2  io  m 

^ 

h 

o 
CD 

E 

13 

X 

Heat   Transfe 

Coefficient   Multipl 

Frequency  =  31. 2 

Diameter  =  3.1* 

% 

< 

^\ 

°o- 

IT) 


I 


CM 


ro 


E 
o 

(/> 

c 


LU 

or 

id 

c/) 

LU 

or 
o_ 
or 

LU 

> 

o 


00 


o 

c 

•r-l 

e 

u 

<D 

Q 

0) 

3 


Oh 

S 
< 
t— i 

o 
•t-< 

Oh 

H 


tM 

u 


o 

cvj 


00 


<s> 


CvJ 


CO 

d 


>W0M    01    3>1    11X3    -OllVd 


the  estimate  is  high, the  ratio  exceeds  unity.     Numerical  interpolation 
helps  determine  the  final  calculated  overpressure. 

3.     General  Theory  of  Frequency  Calculations 

Except  for  temperature  gradient  and  heat  transfer  effects,    the 
pressure  pulses  of  thermal  oscillations  are  the  same  as  those  of  water 
hammer.     As  described  by  Rich  (1951),    one  common  case  of  water  ham- 
mer involves  pressure  pulses  traveling  in  a  pipe  with  its  open  end  in  a 
large  reservoir  and  its  closed  end  at  a  valve.     In  that  case  the  initial 
condition  is  liquid  flowing  in  a  long  pipe  from  a  reservoir  through  a  valve 
to  equipment  of  some  sort  or  to  discharge.     If  the  valve  is  suddenly 
closed,    there  is  a  sequence  of  four  steps  which  cyclically  repeat  them- 
selves until  damped  out  by  friction.      These  steps  are  as  follows. 

Step  one  starts  just  after  the  valve  closure.      Liquid  velocity  in  the 
immediate  vicinity  of  the  closure  is  transformed  into  pressure.      This 
pressure  pulse  travels  to  the  open  reservoir  end  at  sonic  speed.     When  it 
arrives  at  the  reservoir,    the  pipe  is  momentarily  filled  with  stationary 
liquid  whose  pressure  is  higher  than  that  in  the  reservoir.      The  amount 
of  pressure  rise  in  step  one  is  given  by 

AP    =     Cp    |U| 

where 

2 
AP  =    pressure  rise    -------     dynes/cm 

C       =     sonic  velocity    -------     cm/sec 

3 

p        =     fluid  density  --------     g/cm 

|U|     =     magnitude  of  fluid  velocity-    cm/sec. 

Step  two  starts  just  when  the  overpressure  wave  of  step  one 
reaches  the  reservoir.  The  overpressure  at  the  open  end  drops  to 
reservoir  pressure.     In  the  process  it  acquires  a  velocity  equal  to  its 


original  velocity  but  in  the  opposite  direction.      This  reverse  velocity  and 
reservoir  pressure  wave  travels  from  the  reservoir  to  the  closure  at 
sonic  speed.     When  it  arrives  at  the  closed  end  the  pipe  is  momentarily 
filled  with  liquid  at  reservoir  pressure,   whose  velocity  is  everywhere 
the  reverse  of  the  original. 

Step  three  starts  as  soon  as  the  reverse  velocity  wave  arrives  at 
the  closed  end.      Liquid  velocity  in  the  immediate  vicinity  of  the  closure 
is  transformed  into  pressure  drop.     This  reduced  pressure  pulse  travels 
to  the  open  reservoir  end  at  sonic  speed.     When  it  arrives,    the  pipe  is 
momentarily  filled  with  stationary  liquid  whose  pressure  is  lower  than  that 
in  the  reservoir.      The  amount  of  pressure  drop  is  the  same  as  the  pres- 
sure rise  of  step  one. 

Step  four  starts  just  when  the  underpressure  wave  of  step  three 
reaches  the  reservoir.     The  underpressure  at  the  open  end  rises  to  reser- 
voir pressure  and  acquires  the  original  velocity.      This  pressure  and 
velocity  travels  from  the  reservoir  to  the  closure  at  sonic  speed.     Momen- 
tarily at  the  end  of  step  four  the  liquid  pressure  and  velocity  are  the  same 
as  they  were  at  the  beginning  of  step  one. 

After  step  four  the    process  repeats  itself  until  friction  damps  out 
the  oscillations.. 

Except  for  the  effects  of  heat  transfer  and  gas  temperature  gradient, 
the   process  is  similar  for  thermal  oscillations.     Because  of  liquid  evapor- 
ation or  other  causes,    a  small  pressure  rise  occurs  in  the  dewar.     Momen- 
tarily the  tube  is  filled  with  stationary  gas  whose  pressure  is  lower  than 
that  in  the  dewar.      This   is  the  same  situation  that  exists  in  the  pipe  in  the 
water  hammer  case  at  the  end  oi  step  three.     Again  ignoring  the  effects  of 
heat  transfer  and  gas  temperature  gradient,    the  thermal  oscillations 
follow  steps  four,    one,    two,    three  of  the  water  hammer  case. 


If  thermal  oscillations  followed  the  foregoing  sequence  of  events 
exactly,    a  comparatively  simple  theory  ot  heat  transfer  from  the  warm 
regions  of  the  tube  to  the  dewar  could  be  derived.     Even  though  the  actual 
case  is  more  complicated,    an  examination  of  the  simple  case  is  instructive, 

During  step  four  every  particle  of  gas  in  the  tube  is  pressurized  a 
small  amount  and  travels  a  short  distance  into  a  warmer  region  of  the 
tube.     Some  heat  is  transferred  after  pressurization.     During  step  one 
the  gas  is  pressurized  once  more     and  again  moves  into  a  still  warmer 
region.     More  heat  is  transferred.     During  step  two  the   gas  is  depres- 
surized  and  moves  into  a  cooler  region.     Some  heat  is  transferred.     Once 
equilibrium  has  been  obtained  the  final  condition  at  the  end  of  step  three 
and  the  original  condition  at  the  start  of  step  four  are  identical. 

This  chain  of  events  is  shown  on  temperature,    entropy  coordinates 
on  figure  1 . 

A  fluid  which  follows  the  path  shown  must  do  work  on  the  surround- 
ings.    Part  of  the  work  done  is  that  required  to  overcome  wall  and  internal 
gas  friction  to  maintain  the  process.      The  balance  oi  the  work  done  is  ex- 
pended in  expelling  gas  from  the  tube  to  the  dewar.      This   gas  expelled  into 
the  dewar  contains  the  energy  that  increases  the  liquid  evaporation  rate. 

If  the  behavior  oi  a  thermal  oscillation  were  exactly  the  same  as 
that  of  water  hammer,    the  frequency  could  be  determined  by  dividing  the 
tube  into  regions  of  known  temperature,    determining  the  average  sonic 
velocity  in  each  region,    and  calculating  the  time  required  for  a  sonic 
pulse  to  travel  from  one  end  of  the  tube  to  the  other.      The  period  of  the 
oscillation  by  this  method  would  be  the  time  required  for  the  pulse  to 
traverse  the  tube  four  times,    and  the  frequency  would  be  the  reciprocal 
of  the  period.      This  method  was  applied  to  the  data  of  Bannister  (1965). 
As  is   shown  in  table  I,    the  method  is  not  sufficiently  accurate. 


Position 

Length 

Temp. 

Sonic  Vel 

Inches 

cm 

°K 

cm/sec 

0-8 

20.3 

4.  3 

12.  2     x  10^ 

14.4  x  10 
20.  38  x  10 
55.8     x  10 

91. 5  x  10 
101.  0     x  10 

8-16 

2  0.  3 

6 

16-24 

20.  3 

12 

24  -  32 

20.  3 

90 

32  -  40 

2  0.  3 

242 

40  -  48 

20.  3 

295 

TABLE  I 

Using  Temperatures  of  Figure  2 

Time 
sec 

1.664  x  10";? 

1.411  x  10" 

.999  x  10" 

. 364  x  10" 

. 222  x  10" 

.201  x  10" 

-3 
Time  for  one  traverse  4.  861  x  10      sec 

Period  0.  01944  sec 

Frequency  calculated  51.4cps 

Frequency  measured  34  cps 


Since  the  calculated  frequency  of  table  I  does  not  agree  with  meas- 
urement,   it  is  evident  that  the  method  is  inadequate.     The  reason  is 
simply  that  a  temperature  gradient  exists  in  the  pipe  and  causes  pressure 
wave  attenuations  and  reflections.      Therefore,    the  peak  pressure  at  any 
point  is  a  summation  of  multitudinous   reflections      and  does  not  necessarily 
advance  at  sonic  velocity,    even  though  the  first  warning  of  a  pressure 
pulse  does  travel  at  sonic  speed.      To  calculate  frequency  it  is  found  con- 
venient to  divide  the  pipe  into  a  large  number  of  constant  temperature 
regions,    assume  the  initiation     of  a  small  pulse  in  the  dewar,and  calculate 
the  resulting  pressure  versus  time  throughout  the  pipe  from  equations 
derived  in  Appendix  A. 

In  the  derivations,    all  pressures  are  considered  as  absolute.     How- 
ever,   since  the  derived  expressions  show  pressure  differences,    it  is  often 
possible  in  the  digital  computer  program  to  ignore  the   absolute  pressure 
level  and  to  consider  only  over  or  underpressure.      Of  course,    in  deter- 
mining gas  densities,    the  absolute  pressure  must  be  used. 


10 


These  derived  relationships  have  been  programmed  for  calculation 
on  a  high  speed  digital  computer.     As  a  first  operation  the  program  divides 
the  pipe  into  gas  elements  such  that  each  has  an  equal  sonic  traverse  time. 
The  number  of  these  elements  may  be  varied.      The  length  of  time  required 
for  sound  to  traverse  one  tube  length  is   one  calculation  time.     Temperature 
gradients  are  obtained  from  Bannister  (1965).     The  original  pressure  level 

is  an  input  variable.      The  amplitude  of  the  assumed  pressure  pulse  is 

2 

arbitrarily  assumed  to  be  1000  dynes/cm    ,    although  other  values  may  be 

used.     For  as  many  calculation  times  as  desired,    the  program  calculates 
and  prints  gas  velocities  and  pressures  for  each  of  the  N  elements.     By 
examination  of  the  printout    it  is  possible  to  determine  the  pulse  frequency. 
The  program  presented  herein  determines  the  period  between  each  of  the 
first  10  pressure  minimums  thereby  obtaining  9  periods,    and  9  correspond- 
ing frequencies.      The  average  of  these  9  frequencies  is  then  considered 
to  be  the  true  frequency.     Minor  modifications  can  increase  the  number 
of  individual  periods  and  frequencies  averaged. 

4„     Description  of    Amplitude  Calculations 

For  amplitude  calculations   some  assumptions  were  made  as  to  the 
nature  of  the  pressure  oscillations.     It  was  noticed  during  frequency  cal- 
culations that  even  a  step  pressure  input  resulted  in  almost  sinusoidal 
pressure  variation  at  the  closed  end.     Furthermore,    all  peak  pressures 
tended  to  occur  at  the  same  time,    although  the  amplitude  varied  from  end 
to  end.      From  oscilloscope  observations   it  has  been  noticed  that  the  oscil- 
lations tended  towards   sinusoidal.      Therefore,    it  is  assumed  that  at  any 
position  in  the  pipe  the  pressure  oscillation  is  sinusoidal.     The  end  to  end 
variation  of  the  maximum  overpressure  is  obtained  from  a  simple  cosine 
curve  and  the  fact  that  at  the  open  end  there  is  a  pressure  node  at  which 
the  maximum  overpressure  is  always  zero. 
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For  amplitude  determinations  the  duct  is  divided  into  a  large 
number  of  elements,    all  of  which  have  the   same  length.      To  account  for 
gas  entering  and  leaving  the  pipe  an  additional  quantity  of  elements  is 
assumed  to  exist  in  the  dewar.      Thus,    there  are  N  elements  in  the  pipe 
and  M  elements  altogether.      The  period  of  one  cycle  is  divided  into  360 
degrees.      One  time  interval  is  taken  to  be  a  small  fraction  of  the  total. 
Usually,    a  5  degree  time  interval  is   selected  as  being   small  enough  for 
accuracy  without  requiring  excessive  computer  cost. 

Because  of  pressure  change  in  the  element  nearest  the  closed  end, 
element  one,    its  midpoint  will  oscillate  back  and  forth  in  the  duct.     Its 
temperature  will  change  due  to  compression  and  heat  transfer.      The  point 
on  the  wall  in  contact  with  the  midpoint  will  change  with  the  movement, 
resulting  in  temperature  difference  and  heat  transfer.      There  also  will  be 
density  and  position  changes  that  are   somewhat  different  from  an  adiabatic 
case. 

Initial  approximate  temperatures  and  densities  are  assumed. 
After  heat  transfer  has  been  calculated  repeatedly  around  one  complete 
cycle,    conditions  are  different  than  initially  assumed.      However,    after 
a   sufficient  number  of  calculation  cycles,    equilibrium  is  obtained.      Then 
the  position  of  the  cold  end  of  element  one  is  a   series  of  points  varying 
in  a  cycle,    but  invariant  from  one  cycle  to  the  next. 

This  final  series  of  positions  of  the  cold  side  of  element  one  is  the 
permanent  series  of  positions  of  the  warm  side  of  element  two.      Initial 
approximate  temperatures  and  densities  are  assumed  in  element  two  fol- 
lowed by  the   same  type  of  heat  transfer  calculations  that  were  performed 
for  element  one.     After  numerous  cycles  of  calculations,    a  series  of  posi- 
tions is  obtained  for  the  cold  side  of  element  two.      These  are  the  permanent 
series  of  positions  of  the  warm  side  of  element  three. 
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In  general  this  same  series  of  calculations  is  performed  for  each 
of  the  M  elements.     In  each  time  interval  the  heat  transfer,    DQ,    and  the 
infinitesimal  work  done,    DWK,    are  computed.      The  summation  of  these 
quantities      during  the  final  cycle      gives  the  net  heat  transfer  into,and 
the  net  work  done  by, the  I'th  element,    Q(I)  and  WORK(I).     Because  of 
the  oscillation,    gas  is  continually  entering  and  leaving  the  duct.     As  a  re- 
sult there  are  three  types  of  elements.     One  type,    like  element  number  one, 
is  entirely  inside  the  duct  both  at  the  beginning  and  at  the  end  of  a  time 
interval.     A  second  type  is  partly  inside  the  duct  and  partly  in  the  dewar, 
either  at  the  beginning  or  at  the  end  of  the  time  interval,    or  both.     The 
third  type  is  out  of  the  duct  and  inside  the  dewar  at  all  times. 

For  the  first  type,   the  calculations,    though  involved,    are  relatively 
straightforward.       For  the  third  type,    practically  no  calculations  are 
needed.     For  the  second  type  it  is  necessary  to  divide  the  gas  element  into 
a  section  inside  the  pipe  and  a  section  outside  the  pipe,    at  the  beginning 
and  at  the  end  of  the  time  interval.     As  a  result  there  are  three  routes 
through  the  amplitude  program,    one  for  each  type  of  calculation. 

Using  Fortran  notation  the  cosine  variation  of  overpressure  with 
distance  along  the  pipe  is  expressed  as, 


where 


PX(I)  =     PM*COSF(PI/2.    *(XAVG(2)-XX(1)/XL) 

2 
PM  =     Maximum  overpressure  at  closed  end  -  dynes/cm 

PI  =    3.1415926536 

XAVG(Z)      =     Position  of  element  center   -  cm 

XL  =     Length  of  duct  -  cm 

PX(I)  =     Maximum  overpressure  in  element  (I)  at 

2 
position  XAVG(2)  -  dynes/cm 
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XX(1)    =     Position  of  closed  end  of  duct   -  cm. 
The  above  relationship  holds  if  (XAVG(2)-XX(  1))  is  less  than  XL.     Other- 
wise PX(I)  becomes  zero. 

The  momentary  overpressure  at  the  end  of  a  given  time  interval 
in  the  cycle  is 

P(2)  =  PX(I)*SIN2 

where 

2 
P(2  )     =     overpressure  at  end  of  interval  -  dynes/cm 

SIN2    =     the  sine  function  of  the  fraction  of  the  cycle  completed 

at  the  end  of  the  interval. 

One  of  the  primary  variables  calculated  during  each  time  interval 
is  the  infinitesimal  work  done.      The  obvious  way  to  make  this  calculation 
is  to  multiply  the  average  pressure  during  the  interval  by  the  volume  in- 
crease during  the  interval.      In  Fortran  notation  this  is 


DWK    =     PRESS*(VOL(2)-VOL(l)) 


where 


DWK       =     infinitesimal  work  done  in  the  interval  -  ergs 

PRESS    =     average  absolute  pressure  in  the  interval  - 

2 
dynes/cm 

3 

VOL(2)  =     element  volume  at  end  of  the  interval  -  cm 

3 
VOL(l)  =     element  volume  at  beginning  of  the  interval  -  cm    . 

An  instability  in  the  program  was  finally  traced  to  this  method  of  com- 
puting work.     An  equivalent  expression,    derived  in  Appendix  B,    was 
finally  used.      The  equivalent  method  produced  infinitesimal  work  terms 
that  at  equilibrium  were  the  same  as  those  of  the  original  method  with  an 
error  less  than  one  part  in  5000,    and    the  instability  was  largely  eliminated. 
This  method  of  calculation,    in  Fortran  notation,    is 
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DWK  =  (XK-1.  )/XK*DQ-VOLM*(P(2)-P(l))/XK-DFR/XK 
where 

XK  =     specific  heat  ratio 

DQ  =     infinitesimal  quantity  of  heat  transferred  in  a  time 

interval  -  ergs 

3 
VOLM    =     average  element  volume  during  interval  -  cm 

P(2)  =     element  overpressure  at  end  of  time  interval  - 

dynes/cm 
P(l)         =     element  overpressure  at  beginning  of  interval  - 

dynes/cm 
DFR        =     infinitesimal  frictional  energy  loss  during  interval  - 

ergs. 

Kinetic  energy  is  calculated  in  two  ways.      One  method  is  to  deter 
mine  the  infinitesimal  kinetic  energy  in  a  given  element  every  time  the 
velocity  is  towards  the  dewar,    that  is  to  say  when  it  is  positive.      This  is 
done  with  the  equation 

DKE  =  0.  5*(VEL*AF*ROE*DT)*VEL*VEL 

"where 

DKE     =     infinitesimal  kinetic  energy  -  ergs 

VEL    =     gas  velocity  -  cm/sec 

2 
AF       =     pipe  cross  section  area  -  cm 

3 
ROE     =     average  density   -  g/cm 

DT        =    time  interval  -  sec. 

The  terms  in  the  parentheses  are  the  mass   of  gas  moving  past  a  point 
during  the  interval,    so  the  equation  is  the  usual  one  in  which  kinetic 
energy  is   one  half  the  product  of  the  mass   and  the  velocity  squared. 

By  summing  up  the  DKE  for  each  time  interval,    the  dewar  bound 
kinetic  energy  in  one  cycle  for  each  element,    ZKE(I),    is  found.      The 
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kinetic  energy  of  the  N'th  or  M'th  elements,    ZKE(N)  or  ZKE(M),    should 
be  representative  of  the  kinetic  energy  expelled  into  the  dewar. 

A  more  correct  value  of  the  kinetic  energy  expelled  into  the  dewar 
is  that  at  the  exit  from  the  duct.      During  one  cycle,    during  each  time  in- 
terval,   some  one  element  is  partly  in  and  partly  out  of  the  duct.     By  sum- 
ming the  infinitesimal  amount  of  kinetic  energy  of  these  elements,    a 
close  value  of  the  kinetic  energy  expelled  into  the  tank  is  obtained.      The 
equation  used  is 

EKE(J)  =  0.  5*(VEL*AF*ROE*DT)*VEL*VEL 
where 

EKE(J)    =     infinitesimal  kinetic  energy  at  the  duct  exit  during 
the  J'th  time  interval  -  ergs. 

Kinetic  energy  is  obtained  in  the  program  by  both  methods. 
Approximately  the  same  summations  are  obtained  from  each  method. 
The  values  obtained  from  summation  of  EKE(J)  are  used  because  they 
more  closely  reflect  actual  operation,, 

The  program  for  calculating  amplitude  is  called  THM  LOS.     Since 
this  program  does  nothing  except  call  the  subroutines,    READ  DATA, 
ELEMENT,    and  AMPLITUD,    there  is  no  need  for  discussion  oi  it,    but 
there  follows  a  detail  description  of  the  subroutines. 

5.     Detail  of  Subroutine  READ  DATA  in  Amplitude  Calculations 

This  subroutine  is  used  to  read  information  into  the  program  and 
to  calculate  certain  derived  quantities.      Most  of  the  input  information  is 
printed  out  to  show  the  calculation  conditions.     In  the  following  descrip- 
tion    no  explanation  is  made  of  Fortran  statements  which  merely  have 
to  do  with  the  mechanics  of  the  digital  computer  program     such  as 
COMMON,    TYPE  REAL,    or  FORMAT  statements.     Similarly  no  explanation 
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is  made  of  statements  which  themselves  are  merely  explanations.     Sub- 
routine READ  DATA  is  given  in  Appendix  E. 

In  line  5  a  number  of  quantities  are  read  which  have  the  following 
meanings: 

PSTART  =     dewar  pressure  -  dynes/ cm 

DIAM  =     the  duct  diameter.   (This  is  read  in  inches  and 

converted  to  centimeters) 
FRF  =     Fanning  friction  factor 

FACTOR  =     a  multiplier  used  if  desired  to  alter  text  book 

heat  transfer  coefficients 
AVGFREQ       =     the  frequency  of  the  oscillations  -  cycles/sec 
PM  =     the  assumed  maximum  overpressure  -  dynes/cm 

N  =     the  number  of  elements  in  the  duct 

NJJ  =     the  maximum  number  of  iterations  permissible 

NFR  =     a  multiplier  used  if  desired  to  alter  text  book 

friction  factors 
JD  =     the  number  of  degrees  in  one  time  interval.  (  JD 

can  have  the  values  1,  2,  3,4,  5,  6,  8,  9,  10.     The 
value  5  has  usually  been  used.     One  cycle  is   360 
degrees,    so  JD  is  a  fraction  of  a  cycle). 

Since  the  program  determines  the  oscillations  that  result  from  a 
given  temperature  gradient,   temperature  versus  duct  positions  are  needed. 
This  can  be  from  test  data,   or  from  an  arbitrary  gradient  to  determine  its 
oscillation  characteristics.     In  line  10  the  temperature,    TMP,   in  degrees 
Kelvin  is  read  for  each  of  32  duct  positions,    S,  which  are  read  in  inches. 

Viscosity  of  the  gas,   VIS,   in  micropoise  is  read  at  16  temperatures, 
TMRVIS,   in  degrees  Kelvin  at  line  15. 
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Thermal  conductivity,    CONGAS,    in  milliwatts/cm- °K  is  read 
at  12  temper atures,    TCONGAS,    in  degrees  Kelvin  at  line  20. 

Viscosity  and  conductivity  are  obtained  from  Johnson  (I960)  either 
directly  or  by  extrapolation. 

To  determine  natural  convective  heat  transfer  coefficients,    a 
table  of  Nusselt  number  versus  the  Grashof-Prandtl  product  is  needed. 
The  log        of  the  Nusselt  number,    ZNU,    is  read  at  16  values  of  the  log 
of  the  Grashof-Prandtl  product,    GRPR,    in  line  2  5.     Figure  7-10  from 
McAdams  (1954)  supplied  13  of  these  values.      The  two  largest  values 
used,    not  obtainable  from  figure  7-10,    were  obtained  from  figure  7-7  of 
the  same  reference.     Since  the  Nusselt  number  goes  to  zero  if  the 
Grashof-Prandtl  product  goes  to  zero,    it  was  desired  to  approximate  this 
condition.      To  do  so  the  log        of  the  Nusselt  number  is  set  to  -37  when 
the  log        of  the  Grashof-Prandtl  product  is   -37. 

In  lines   30  and  40012  the  duct  positions,    S,    and  the  duct  diameter, 
DIAM,    are  converted  from  inches  to  centimeters. 

In  lines  35  and  40011  the  viscosity,    VIS,    is  converted  from  micro- 
poise  to  poise. 

In  lines  40  and  40010  the  conductivity,    CONGAS,    is  converted  from 
mw/cm-°Kto  erg/sec -cm- °K. 

In  line  45  the  molecular  weight  of  helium,    WT,    the  universal  gas 
constant,   RU,    and  the  value  for  PI  are  entered.     From  these  the  helium 
gas  constant,    R,    and  the  cross  section  area  of  the  duct,    AF,    are  calcu- 
lated.    In  line   50  the  specific  heat  ratio,    XK,    is  entered,    and  the  specific 
heats,    CPG     and  CVG,    are  calculated.     In  line   55  the  temperature  oi  the 
gas  in  the  dewar  just  outside  the  duct,    TDEW,    is  entered  and  the  density, 
RHODEW,    is  calculated.      The  quantity,    M,    is  also  calculated.     Without 
oscillations  there  are  N  gas   elements  in  the  duct.      With  oscillations, 


elements  or  portions  of  elements  enter  the  duct  from  the  dewar  and 
enter  the  dewar  from  the  duct.      To  keep  track  of  all  these  elements  an 
additional  quantity  of  elements  is  considered.      Thus,    M  represents  a 
number  of  elements  some  of  which  are  always  in  the  duct,    some  are 
partly  in  and  partly  out  of  the  duct,    and  some  are  always  in  the  dewar. 

The  remainder  of  the  subroutine  READ  DATA  is  used  for  printing 
out  the  input  information.     Computational  units  are  printed  except  in  the 
case  of  viscosity,    VIS,   where  convenience  and  space  cause  it  to  be  printed 
in  micropoise  rather  than  poise. 

The  input  value  of  the  Fanning  friction  factor,    FRF,    is  not  printed 
in  READ  DATA.     As  explained  later,    the  input  value  is  subject  to  change. 
The  range  of  values  actually  used  is  printed  from  subroutine  AMPLITUD. 

6.     Detail  of  Subroutine  ELEMENT  in  Amplitude  Calculations 

The  subroutine  ELEMENT,    when  used  with  amplitude  calculations, 
divides  the  duct  into  N  elements  and  M  elements  total  in  the  duct  and 
dewar.     All  elements  are  of  the  same  length,    DXX.     The  position  of  each 
boundary,    XX,    the  wall  temperature  at  each  boundary,    TEMP,    the  initial 
average  temperature  of  each  element,    TMG,    the  initial  average  density  oi 
each  element,   RHO,   the  mass  of  gas  in  each  element,    GM,    the  initial  mid- 
point of  each  element,    XMID,    are  also  calculated  in  ELEMENT  and 
printed.      These  are  all  initial  or  non-oscillating  conditions.     This  sub- 
routine,   ELEMENT,    is  given  in  Appendix  F. 

In  line  5  the  warm  boundary  temperature  of  the  first  element, 
TEMP(l),  is  set  equal  to  the  temperature  read  in  at  that  point,    TMP(l). 
Similarly  the  warm  boundary,    XX(1),    is  set  equal  to  S(l).     The  quantities 
NP1  and  MP1,    needed  later,    are  set  equal    to  N+l  and  M+l.     Line  15  sets 
the  cold  end  boundary,    XX(I+1),    of  elements   1  through  N.     Lines  20  through 
33  compute  the  temperature,    TEMP(I+1),    corresponding  to  each  position, 
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XX(I+1).     Lines   35  and  40  calculate  the  average  initial  element  tempera- 
ture,   TMG(I),    the  average  initial  density,    RHO(I),    the  initial  average 
position,    XMID(I),    and  the  mass  of  gas  in  each  element,    GM(I)„     Lines 
45  through  60  compute  the  same  quantities  for  elements  N+l  through  M. 
The  remainder  of  the  program  prints  out  the  values  of  these  quantities. 

7.     Detail  of  Subroutine  HEAT  TR 

When  called  by  subroutine  AMPLITUD,    subroutine  HEAT  TR  cal- 
culates the  friction  factor,    FRF,   the  Reynolds  number,    REYNO,    the  heat 
transfer  coefficient  by  forced  convection,    HHT1,    the  Grashof-Prandtl 
product,    GP,    the  heat  transfer  coefficient  by  natural  convection,    HHT2, 
the  total  heat  transfer  coefficient,    HHT,    and  the  heat  transferred,    dQ. 
When  it  is  desired  not  to  calculate  the  friction  factor  but  to  use  a  constant 
input  friction  factor,    the  subroutine  omits  calculation  of  friction  factors. 
Subroutine  HEAT  TR  is  given  in  Appendix  G. 

In  line  1,  if  the  variable,  ITESTFRF,  has  been  set  equal  to  unity 
in  subroutine  AMPLITUD,  only  friction  factors  are  desired  and  the  pro- 
gram jumps  to  line  386. 

Lines  2  through  3  calculate  the  momentary  conductivity,    CONG,    of 
a  gas  element  from  the  momentary  element  temperature,    TAVG,    by 
linear  interpolation  in  tables  which  have  been  read  in  the  program. 

Lines  4  through  5  calculate  the  momentary  viscosity,    VISC,    of  the 
gas  element  from  the  momentary  element  temperature,    TAVG,   also  by 
linear  interpolation. 

At  line  6,    if  the  viscosity  is  zero,    an  error  exists  and  the  program 
is  terminated. 

In  line  8  the  Reynolds  number,    REYNO,    is  calculated  using  vari- 
ables determined  in  AMPLITUD.      Since  a  Reynolds  number  is  always 
positive,    its  absolute  value  is  used. 
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In  line  10  the  lowest  possible  value  of  the  Reynolds  number  is  set 
to   0.  001  to  avoid  the  possibility  of  division  by  zero.     Furthermore,    when 
the  Reynolds  number  is  less  than  0.  001  forced  convection  is  usually 
negligible-     Friction  may  be  appreciable,   but  this  assumption  merely 
forces  the  friction  factor  to  a  single  high  value  when  the  Reynolds  number 
would  be  below  0.  001. 

In  line  22,   if  the  conductivity  of  the  gas  is  negative  or  zero, an 
error  exists  and  the  program  terminates. 

In  line  25  the  Prandtl  number  is  calculated.     If  it  is  zero  or  negative, 
an  error  exists  and  line  26  initiates  program  termination. 

In  line  30,   if  the  momentary  gas  temperature,    TAVG,   is  negative 
or  zero, an  error  exists  and  the  program  terminates. 

Line  31  calculates  the  Nusselt  number,    ZNUSSELT,  by  equation 
(9-10a)  of  McAdams(1954).     From  this,    line  32  determines  the  forced 
convective  heat  transfer  factor,    HHT1.     The  multiplier,    FACTOR,   is  an 
input  quantity  used  to  determine  the  effect  on  oscillations  of  reduced  or 
increased  heat  transfer  factors. 

Line  33  determines  the  heat  driving  temperature  difference, 
THETA.     If  it  is  zero,    the  program  skips  unnecessary  calculations  and 
goes  to  line  385. 

Lines  35  through  42  calculate  the  log         of  the  momentary  Grashof- 
Prandtl  product,   GPL.     If  the  Grashof-Prandtl  product,    GP,    is  negative, 
an  error  has  occurred  and  the  program  is  sent  to  line  410  and  terminated. 
If  it  is  zero  the  program  goes  to  line  380  and  skips  unnecessary  calculations. 

Lines  43  through  370  determine  the  login    °"^  the  free  convective 
Nusselt  number,    ZNUL,    and  line  375  determines  the  corresponding 
Nusselt  number,    ZNUS. 
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In  line  380  the  free  convective  heat  transfer  coefficient,    HHT2,    is 
determined.      The  multiplier,    FACTOR,    is  an  input  quantity  used  to  deter- 
mine the  effect  of  increased  or  decreased  heat  transfer  factors  on  oscil- 
lations.     The    overall  heat  transfer  coefficient,    HHT,    is  obtained  by- 
simple  addition  of  forced  and  natural  heat  transfer  coefficients. 

In  line  385  the  heat  transfer  coefficient,    HHT,    is  multiplied 
by  the  heat  transfer  area,    AREA,    the  temperature  difference,    THETA, 
and  the  time  increment,    DT,    to  obtain  the  infinitesimal  heat  transferred, 
DQ,    for  the  element  in  the  time  increment. 

Lines   386  through  39  5  have  to  do  with  friction  factor.     In  line  386 
the  minimum  value  of  the  Reynolds  number  is  limited  to  0.  001. 

In  line  2386  the  decision  is  made  whether  or  not  to  calculate  the 
friction  factor.      The  program  has  the  option  of  using  a  constant  input 
friction  factor  or  oi  using  a  modified  text  book  friction  factor.     If  the 
value  of  the  friction  factor  during  input  is  zero,    then  in  subroutine 
AMPLITUD  the  variable,  IFRF,    is  set  to  zero.     Alternately,    if  the  input 
value  of  the  friction  factor  is  not  zero,    the  variable,    IFRF,    is  set  to  unity. 
Thus,    line  2386  causes  calculation  of  a  friction  factor  if  the  original  input 
friction  factor  is  zero.     Otherwise,    it  sends  the  program  to  line  99999 
and  makes  no  friction  factor  calculations. 

If  the  Reynolds  number  is  less  than  1000  in  line  387,    the  friction 
factor  is  obtained  from  the  usual  laminar  equation,    but  multiplied  by  NFR. 
The  quantity,    NFR,    is  an  input  multiplier  used  to  explore  the  effect  of 
increased  or  reduced  friction.     If  the  Reynolds  number  is  less  than  1000 
in  line  388,    the  friction  factor  has  already  been  obtained,    and  control 
returns  to  AMPLITUD. 

If  the  Reynolds  number  is   greater  than  or  equal  to  1000,    lines 
389  through  39  5  calculate  the  friction  factor  by  an  equation  (6-8e)  of 
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McAdams  (1954).      The  multiplier  NFR  is  used  as  in  laminar  flow.      The 
method  is  iterative.     A  value  for  friction  factor  is  assumed  in  lines   389 
and  390.     A  better  estimate  is  found  in  line  392.     If  this  estimate  dis- 
agrees with  the  previous  estimate  by  an  appreciable  amount  the   program 
returns  to  line  390  and  iterates  again.     Otherwise,    the  calculation  of 
friction  factor  is  complete. 

Lines  400  through  420  are  used  in  case  of  error.  In  the  event  of 
an  unstable  calculation,  TAVG,  CONG,  or  GP  can  become  negative.  If 
this  occurs,  there  is  a  print-out  to  aid  in  diagnosis  and  the  quantity,  DT, 
is  set  equal  to  zero,   which  causes  termination  of  the  program. 

8.     Detail  of  Subroutine  AMPLITUD 

This  subroutine  is  the  fundamental  part  of  the  program.     All  the 
other  subroutines  are  for  the  purpose  of  helping  subroutine  AMPLITUD. 
As  an  element  changes  pressure,    its  density  changes,    heat  is  transferred 
from  or  to  the  wall,   work  is  done  by  the  element,    friction  is  experienced, 
and  at  the  open  end,    kinetic  energy  is  lost  from  the  pipe  and  ejected  into 
the  dewar.     From  an  assumed  value  of  the  maximum  overpressure  at  the 
closed  end,    AMPLITUD  computes  the  behavior  of  each  element  in  turn 
from  1  to  M.     Subroutine  AMPLITUD  is  given  in  Appendix  H. 

The  line  before  line   1  determines  whether  or  not  the  program  will 
use  input  or  calculated  friction  factors.     If  the  input  value  of  friction 
factor,    FRF,    is  zero,the   parameter,    IFRF,    is  set  equal  to  zero.     If  the 
input  value  is  different  from  zero,    the  parameter,    IFRF,    is   set  equal  to 
unity.      Then,    if  the  input  friction  factor  is  zero,    IFRF  is  zero,    and  sub- 
routine HEAT  TR  computes  friction  factors  for  each  element  in  each  time 
interval.     Alternately,    if  the  input  friction  factor  is  any  desired  non-zero 
quantity,    IFRF  is  unity,    and  subroutine  HEAT  TR  leaves  the  friction  factor 
unchanged  at  its  input  value. 
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Lines   1  through  40009  cause  the  printing  of  a  representative  table 
of  friction  factors  used  in  the  program  versus  Reynolds  number.     Where 
a  constant  friction  factor  is  used,   this  same  friction  factor  shows  for  each 
Reynolds  number.      Where  a  calculated  friction  factor  is  used,    the  friction 
factor  calculated  for  the  corresponding  Reynolds  number  is  printed.      The 
friction  factor  printed  includes  the  multiplier  NFR  for  the  calculated  case. 
When  the  input  friction  factor  is  used,    NFR  is  inoperative. 

Lines  5  through  25  set  the  values  of  certain  constants  needed  during 
calculations  or  set  initial  values  for  variables  to  insure  proper  operation 
of  the  program.      The  most  important  of  these  is  the  time  interval,   DT, 
calculated  in  line  10.      The  time  of  one  interval  is  the  time  for  one   cycle 
times  the  fraction  of  a  cycle  used.      Since  the  time  for  one  interval  is  the 
reciprocal  of  the  frequency,   AVG  FREQ,    and  since  the  cycle  fraction  is 
JD/360,    in  Fortran  notation,    the  time  interval  is  given  by 

DT  =  JD/360. /AVG  FREQ. 

Line  3  0  is  a  DO  LOOP  ending  in  line  40000  which  causes  the  calcu- 
lation of  all  quantities  needed  for  each  of  the  elements   1  through  M. 

Lines  31  through  55,   for  each  gas  element  in  turn,    set  initial  values 
for  variables  to  insure  proper  operation. 

In  the  DO  LOOP  from  line  60  to  line  40002,    the  calculation  for  each 
element  is  caused  to  repeat  as  many  times  as  needed.      The  maximum  num- 
ber of  repetitions  is  the  quantity  NJJ,    an  input  quantity.      The  value  of  NJJ, 
in  general, is  kept  high  enough  to      insure  that  conditions  calculated  for  the 
final  cycle  are  substantially  the  same  as  for  the  previous  one.      The  quantity, 
IQ,    is  used  to  indicate  the  adequacy  of  this  convergence.      Later  in  the  pro- 
gram,   if  cycle  convergence  is  adequate,    IQ  is  increased  by  1.     If  this 
happens  twice,    the  accuracy  is  sufficient  and  the  program  proceeds  to  the 
next  element. 
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In  line  6l,   the  maximum  number  of  cycles,    JJ  MAX(I),    used  for 
the  gas  element  I  is  determined.     Since  this  quantity  is  printed,    it  is 
possible  to  see  if  the  input  quantity  NJJ  is  adequate. 

In  line  65  the  quantity  QI  is  set  equal  to  Q(I)  and  WORK  I  to  WORK(I) 
for  later  use  in  determining  adequacy  of  convergence.      Later  the  sum  of 
all  infinitesimal  heat  transfer  to  one  gas  element,    Q(I),    and  the  sum  of 
the  work  terms  in  one  gas  element,    WORK(I),is  obtained.     At  that  time 
QI  and  WORK  I  become  the  heat  and  work  summations  for  the  previous 
cycle  and  a  comparison  is  then  possible. 

Other  calculations  in  lines  65  and  70  set  initial  values  to  insure 
proper  program  operation. 

Line  80,    a  DO  LOOP  ending  in  line  40003,    causes  calculations  for 
one  cycle  to  occur  for  the  element  already  selected  in  line  30.     It  sets  the 
variable  J  equal  to  the  portion  of  a  cycle,    expressed  in  degrees,    completed 
at  the  end  of  this  time  interval.      The  first  time  interval  is  zero  to  JD  de- 
grees (If  JD  is  5  then  the  interval  is  0-5  degrees.  )  and  so  on  until  the 
cycle  is  completed  at  3  60  degrees.     Because  of  line  80,  the  input  variable 
JD  must  be  divisible  into  360.     /lso  computed  in  line  80  is  the  sine  function 
of  the  angular  fraction  of  the  cycle  completed  at  the  end  of  the  time  interval. 

The  quantity,  JLJD,  refers  to  the  time  (or  angular  fraction  of  the 
cycle)  at  the  beginning  of  the  interval  just  as  the  quantity,  J, refers  to  it 
at  the  end  of  the  interval.     In  line  85  it  is  calculated.      At  the  beginning  of 
the  cycle  J  =  JD,    JLJD  =  360,    but  at  any  later  interval  when  J  >  JD, 
JLJD  =  J-JD. 

In  lines  90  through  105  numerous  variables  at  the  beginning  of  the 
current  interval  are  set  equal  to  the  same  variable  at  the  end  of  the  prev- 
ious interval.      Thus,   RO(l),   the  gas  density  at  the  beginning  of  the  interval, 
is  set  equal  to  RO(2),   the  density  at  the  end  of  the  previous  interval.    Later, 
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a  value  of  the  density  at  the  end  of  the  current  interval  is  determined, 
which  then  becomes  RO(2).     The  same  is  true  for  the  volume  of  the  ele- 
ment,   VOL,    the  average  position  of  the  element,    XAVG,    the  overpressure 
of  the  element,    P,    the  absolute  pressure  of  the  element,    PRES,    the  mass 
of  the  gas  of  an  element  that  is  in  the  pipe,    G,   and  the  mass  of  an  element 
that  is  in  the  dewar,    GID2  and  GID1. 

Line   110  determines  if  a  gas  element  is  entirely  outside  the    pipe 
and  in  the  dewar  by  means  of  the  variable,    X.     Variable,    X,    is  similar 
to  XX.     XX  is  computed  in  ELEMENT  and  is  the  non-oscillating  or  initial 
boundary  of  the  various  elements.      Thus,    XX(1)  was  set  equal  to  input 
position  S(l)  and  is  the  closed  end  of  the  pipe.     XX(Z)  is  the  position  of 
the  boundary  between  elements  land  2.     XX(I+1)  is  the  boundary  between 
elements  I  and  1+1.     The  final  boundary  in  the  pipe  is  XX(N+1)  which  is 
equal  to  S(32).     Also  in  line  25,    to  save  space,    the  variable,    XF,   was 
set  equal  to  XX(N+1)  and  represents  the  boundary  between  the  pipe  and 
the  dewar.     Since  XX  is  invariant  with  time  it  is  single  subscripted.     X 
varies  with  time  as  well  as  with  the  element  being  considered  and,    there- 
fore,   is  double  subscripted.      Thus  X(l,  J)  refers  to  the  warm  side  of  a 
given  element  at  time  J,    and  X(2,  J)  refers  to  the  cold  side  at  time  J.      To 
conserve  storage,    the  first  subscript  is  always   1  or  2.      This  necessitates 
a  change  when  shifting  from  one  element  to  the  next.     The  value  of  X  for 
the  previous  element  that  was  the  cold  side  becomes  the  value  for  the 
warm  side  o±  the  next  element.      This  shift  is  accomplished  in  line  3  5  as 
the  program  progresses  from  one  element  to  the  next.     For  element  one, 
the  initial  values  o±  X  for  all  time  increments  are  set  equal  to  XX(1).     So 
control  shifts  from  line   110  to  line   690  whenever  the  warm  side  position 
of  the  element  is  either  at  the  end  of  the  pipe  or  protruding  into  the  dewar. 
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If  the  element  is  partly  or  wholly  in  the  pipe  at  the  beginning  or  end 
of  the  interval,   the  calculation  continues  to  line  120. 

In  line  120,   to  insure  proper  operation  of  the  program,   two  heat 
transfer  parameters  are  set  to  zero. 

The  program  now  calculates  heat  transfer  from  the  wall  to  the  gas. 
Lines  240  through  490  do  this  for  elements  wholly  in  the  pipe.      Lines  240 
through  3  00  plus  lines  500  through  2689  do  it  for  elements  partly  in  the 
pipe  and  partly  in  the  dewar. 

The  heat  transfer  calculation  is  necessarily  iterative.     From  the 
previous  calculations  on  other  elements,  the  warm  end  position  of  the  ele- 
ment is  known.      Then  from  an  approximate  initial  volume  of  the  element 
its  approximate  average  position  is  found.     Heat  is  transferred  as  if  this 
were  its  real  average  position.     A  more  accurate  density  and  volume  are 
calculated.     From  this  new  volume  a  new  average  position  is  found  and 
the  heat  transfer  is  recalculated.     The  process  is  repeated  until  the  heat 
transfer  calculation  stabilizes. 

In  line  240,    the  variable,  KOUNT,has  1  added  to  it  in  order  to 
print  out  later  the  number  of  heat  transfer  iterations  needed  for  conver- 
gence.     KOUNT  was  initially  set  to  zero  in  line   105  so  that  its  value  always 
represents  the  number  of  heat  transfer  iterations  for  a  given  element  and 
time  interval  during  the  last  cycle  iteration. 

In  line  240,   the  parameter, DQ1, is  set  equal  to  AQ,    the  heat  trans- 
fer from  the  previous  iteration,   for  convergence  test  purposes.     Note: 
AQ  and  DQ  are  essentially  the  same. 

In  line  241,   the  maximum  value  of  KOUNT  for  the  given  element, 
KOUNT  MAX(I),   is  found.     KOUNT  MAX(I)  is  set  to  zero  in  line  65  so  that 
the  final  KOUNT  MAX(I)  is  for  the  final  cycle. 
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Line  245  sets  the  position  of  the  cold  end  of  the  given  element 
equal  to  the  warm  end  position  plus  the  length  of  the  element.     This  position 
varies  slightly  with  each  KOUNT  because  the  heat  transfer  changes  the 
volume  of  the  element.     That  is,    VOL(2)  is  the  element  volume  at  the  end 
of  the  time  interval.     Dividing  by  the  cross  section  area,   AF,    gives  the 
length.     But  VOL(2)  is  not  known  absolutely  until  after  the  heat  transfer 
and  work  terms  are  evaluated.     Initially,    an  approximation  for  VOL(2)  is 
found  in  line  45.      Successive  calculations  make  it  more  accurate.     During 
the  final  heat  transfer  iteration  of  the  final  cycle  the  cold  side  position  of 
the  given  element  at  time  J,    X(2,  J),    is  firmly  established. 

Line  250  determines  the  average  position  of  the  element  at  the  end 

of  the  time  interval,    J,    by  averaging  the  cold  and  warm  end  positions. 

Line  255  establishes  the  gas  velocity  during  the  interval.     It  divides 
the  motion  of  the  average  position  by  the  time  interval. 

Line  260  determines  the  average  gas  temperature  during  the  inter- 
val.    For  an  initial  approximation,    line  35  sets  the  instantaneous  tempera- 
tures during  a  cycle  to  the  gas  temperature  found  in  ELEMENT  for  the 
given  gas  element.     During  iteration  a  more  accurate  value  is  found  for 
the  temperature  at  the  end  of  the  interval,      T(J),    and,    therefore,    for  the 
average  temperature,    TAVG.     The  temperature  determined  at  the  end  of 
the  previous  interval  establishes  the  temperature  at  the  beginning  of  the 
current  interval,    T(JLJD). 

Lines  265  and  270  calculate  the  maximum  overpressure  occurring 
at  the  end  of  the  time  interval  when  the  average  position  is  at  XAVG(2). 
Line  275  multiplies  this  maximum  overpressure  by  the  sine  of  the  angular 
position  to  obtain  the  instantaneous  overpressure  at  the  end  of  the  interval. 
The  corresponding  absolute  pressure,    PRES(2),    is  obtained  by  adding 
dewar  pressure  to  the  overpressure. 
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Line  280  obtains  the  average  absolute  pressure,    PRESS,    and  the 
average  density,    ROE,    during  the  interval  by  averaging  conditions  at  the 
beginning  and  end  of  the  interval. 

Line  300  transfers  calculations  to  line  500  if  the  element  protrudes 
into  the  dewar  at  the  beginning  or  at  the  end  of  the  interval. 

If  the  element  is  entirely  in  the  pipe,    line  305  calculates  its  average 
volume,    VOLM,   by  dividing  the  element  mass,    GM(I),    by  its  average 
density,   ROE. 

Line  405  calculates  the  average  position,    AVGX,    of  the  element 
during  the  interval     by  taking  the  average  of  the  average  positions  at  the 
beginning  and  end  ofthe  interval. 

Lines  410  through  420  calculate  the  wall  temperature,    TWI,    corres- 
ponding to  the  average  position  of  the  element  during  the  interval. 

Line  425  calls  the  subroutine  HEAT  TR.      The  quantity  ITEST  FRF 
is  set  to  zero  so  that  heat  transfer  will  be  calculated.     Had  it  been  set  to 
unity, only  friction  factor  information  would  have  been  returned. 

Sometimes  the  program  is  unstable  and  troubles  occur  in  subrou- 
tine HEAT  TR.     If  this  occurs,    DT  is  set  to  zero  in  the  subroutine  and 
line  430  causes  termination  of  the  program. 

Line  435  calculates  the  infinitesimal  friction  work,    DFR,    done 
in  the  interval  by  a  formula  derived  in  Appendix  C.     Since  the  friction  is 
always  a  positive  number,   the  absolute  value  is  used. 

Line  440  calculates  the  infinitesimal  net  work,    DWK,    after  sub- 
tracting friction  work,    by  a  formula  derived  in  Appendix  B. 

Line  455  calculates  the  internal  energy  change  of  the  element,    QW, 
from  the  first  law  of  thermodynamics  by  subtracting  the  net  work,    DWK, 
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from  the  heat  input,    DQ,    returned  from  HEAT  TR.      The  quantity  AQ, 
used  later, is  set  equal  to  DQ. 

The  gas  temperature  at  the  end  of  the  time  interval,    T(J),    is  cal- 
culated in  line  460.      To  the  temperature  at  the  beginning,    T(JLJD),   is 
added     the  temperature  rise  computed  by  dividing  the  internal  energy 
increase,    QW,   by  the  mass  of  gas  in  the  element,    GM(I),  and  the  specific 
heat,    CVG. 

Lines  465  and  47  0  calculate  the  element  density  and  volume  at  the 
end  of  the  interval. 

During  the  first  iteration  of  a  heat  transfer  calculation  in  an  inter- 
val,   it  is  not  possible  for  equilibrium  to  have  been  established,    so  line 
475  forces  the  program  back  to  line  240  for  at  least  one  more  iteration. 

Sometimes  on  undue  number  of  heat  transfer  iterations  occur  and  it 
is  necessary  to  proceed  without  complete  convergence.     If  there  have  been 
more  than  10  attempts,    KOUNT  will  exceed  10,    and  line  480  sends  the 
program  to  line  700  to  continue  on  with  the  program.      This  seldom  happens. 
If  it  does,   the  quantity, KOUNT  MAX(I), calculated  in  line  241  exceeds  10 
and  shows  in  the  print-out. 

A  test  for  heat  transfer  convergence  is  in  line  485.     If  the  heat 
transfer  computed  is  appreciably  different  from  the  previous  computation, 
the  program  returns  to  line  240  and  tries  again.      To  avoid  division  by  zero 
the  quantity  AQ,    if  zero,    is  changed  to  0.  000001. 

If  the  program  gets  to  line  490, it  has  completed  heat  transfer  and 
work  calculations  for  this  time  interval  of  this  element, and  the  calculations 
have  been  for  an  element  entirely  in  the  pipe.     It  now  skips  calculations 
designed  for  elements  partly  or  wholly  in  the  dewar  and  control  is  shifted 
to  line  7  00. 
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If  the  element  is  partly  in  the  dewar  and  partly  in  the  pipe,    either  at 
the  beginning  or  the  end  of  the  interval,   control  is  transferred  to  line    500  by 

line  300.  For  heat  transfer  purposes  it  is  necessary  to  determine  the  length 
of  the  element  touching  the  wall.     Four  quantities  are  defined,  CI,    C2,    Wl, 
and  W2.      These  are  the  positions  of  the  cold  side  at  the  interval  beginning, 
the  cold  side  position  at  the  interval  end,    and  the  warm  side  at  the  begin- 
ning and  end.     Depending  on  when  and  by  how  much  the  element  protrudes 
into  the  dewar,    these  quantities  may  be  the  sides  of  the  element  or  the 
ends  of  the  pipe.      Lines  500  through  530  determine  these  positions. 

Lines  550  and  551  calculate  the  infinitesimal  amount  of  kinetic 
energy  ejected  into  the  dewar  during  this  time  interval. 

Line  600  calculates  the  heat  transfer  area  at  the  start,    AREA1, 
at  the  end,   AREA2,    and  the  average  area,    AREA. 

Line  605  calculates  the  average  position  of  the  element  for  heat 
transfer  purposes. 

Lines  610  through  620  determine  the  wall  temperature,    TWI,    at 
the    average  element  position. 

Line  625  determines  heat  transfer  into  the  element  by  calling  sub- 
routine HEAT  TR.      The  quantity  ITEST  FRF  is  set  equal  to  zero  so  that 
HEAT  TR  returns  heat  transfer  as  well  as  friction  factor  information. 

Line  63  0  terminates  the  program  if  a  calculation  instability  is 
found  in  HEAT  TR. 

Line  635  calculates  the  infinitesimal  friction  work  done  during  the 
interval  by  a  formula  derived  in  Appendix  C. 

Line  639  calculates  the  overage  volume  of  the  element  in  the  pipe. 
Since  it  is  the  volume  in  the  pipe  rather  than  the  entire  element  volume  that 
is  desired,    the  quantities  G(l)  and  G(2)  are  used  for  the  mass  of  gas  instead 
of  GM(I). 

31 


Line  640  calculates  the  net  work  done  in  the  interval  by  a  formula 
derived  in  Appendix  B. 

Line  655  calculates  the  internal  energy  change  in  the  element. 

Line  670  calculates  the  gas  volume  of  the  given  element  that  is  in 
the  pipe,   VOLIPZ,    the  mass  of  this  gas,    G(2),    and  the  mass  of  gas  of  this 
element  in  the  dewar,    GID2,    at  the  end  of  the  interval. 

Lines  674  and  67  5  calculate  the  density  of  the  fraction  of  the  gas 
element     that  is  in  the  dewar  at  the  end  of  the  interval,    ROID2.      This 
quantity  is  used  in  line  67  6  to  determine  the  volume  of  the  dewar  fraction 
of  the  gas  element,    VOLID2,    and  the  total  volume  of  the  element,    VOL(2). 
Note:     The  calculations  oflines  674through  676  possibly  are  not  completely 
rigorous,   but  they  do  not  affect  the  results  obtained  using  kinetic  energies, 
EKE. 

Lines  68  0  through  687  calculate  the  temperature  of  that  portion  of 
the  gas  element  that  is  in  the  pipe  at  the  interval  end.      Three  routines  are 
necessary,    one  for  flow  out  of  the  pipe,    one  for  flow  into  it,    and  one  special 
case  where  there  is  no  gas  in  the  pipe  at  the  end  of  the  interval.     It  should 
be  noted  that  T(J)  is  always  the  temperature  at  time  J  of  the  gas  fraction 
in  the  pipe. 

Line  688  computes  the  density  of  the  gas  fraction  in  the  pipe  at  the 
interval  end. 

In  lines  688  through  689  criteria  similar  to  those  of  lines  475 
through  485  are  used  to  test  heat  transfer  convergence. 

At  line  2689,   the  program  skips  calculations  used  for  elements  that 
are  entirely  in  the  dewar. 
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Lines  690  through  697  are  for  the  case  of  the  gas  element  entirely 
in  the  dewar.     No  computations  are  really  needed.     Those  that  are  made 
are  for  the  purpose  of  maintaining  consistency  of  printout. 

At  line  7  00  all  calculations  have  been  performed  on  the  gas  element 
for  the  current  time  interval.     In  line  701,    the  infinitesimal  work,    DWK, 
just  calculated,is  added  to  the  work  term  for  the  element  so  that  at  the 
end  of  the  cycle  the  term,  WORK(I),  will  represent  the  total  work  done  by 
the  element  in  the  interval.     Similar  steps  are  taken  to  obtain  the  total 
heat  input  to  the  element,    Q(I),    and  the  total  friction  in  the  element,    FR(I}. 
In  lines  703  and  7  04  the  dewar  bound  kinetic  energy  of  the  element,    ZKE(I), 
is  computed. 

If  the  convergence  from  cycle  to  cycle  is  not  yet  adequate,    the 
quantity  IQ  will  be  zero  and  detail  printing  is  not  desired  so  line  7  05  sends 
the  program  to  line  40003  for  return  to  line  80  to  advance  to  the  next  time 
interval,   or  if  a  cycle  has  been  completed,   to  line  723. 

If  only  one  cycle  has  been  computed  for  the  element,    line  723 
transfers  control  to  line  40002  which  returns  control  to  line  60  to 
start  a  second  cycle. 

Lines  724  through  750  concern  cycle  convergence,    printing,    and 
program  termination. 

There  are  three  criteria  for  cycle  convergence.     For  the  elements 
that  are  wholly  in  the  pipe,   the  heat  received  by  the  element  must  equal 
the  net  work  done  by  the  element.     For  all  elements     at  equilibrium,   the 
heat  received  by  an  element  in  a  calculation  cycle  must  equal  the  heat  re- 
ceived in  the  previous  calculation  cycle.      The  work  done  in  a  calculation 
cycle  must  equal  that  done  in  the  previous  one.      These  three  criteria  are 
calculated  in  the  form  of  three  ratios,   QLWOQ,    DELQIOQI,    and  DELWKOWK. 
They  are  defined  as  follows.     QLWOQ  is  the  difference  between  the  heat  and 
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work,    divided  by  the  heat  input.  DELQIOQI  is  the  difference  between  this 
cycle's  and  the  previous  cycle's  heat,    divided  by  this  cycle's  heat  transfer. 
DELWKOWK  is  the  same  applied  to  work  done.     In  all  three  criteria  the 
absolute  value  is  used. 

If  the  heat  transfer  and  the  work  done  are  both  zero, the  criterion 
is   satisfied  and  line  725  sends  the  program  to  line  750. 

Line  726  avoids  division  by  zero. 

Line  730  and  73  5  calculate  the  three  criteria. 

If  the  values  of  any  two  of  the  three  criteria     are  less  than  0.  0001 
in  lines  738  through  740,    the  program  is   sent  to  line  750. 

The  quantity  NJJL2  was  calculated  in  line  6  and  is  equal  to  two  less 
than  the  input  quantity  NJJ,    which  limits  the  number  of  cycles  that  can  be 
calculated.     If  the  number  of  cycles  is  as  great  as  NJJL2,  it  will  be  desired 
to  print  detailed  information  before  proceeding  to  the  next  element.     The 
program  is  therefore  sent  to  line  750  by  line  742. 

Line  7  50  adds   1  to  the  variable,    IQ.     Then  during  the  next  cycle 
various  details  will  be  printed.     If  cycle  convergency  has  occured  twice, 
the  variable  IQ  increases  to  two  and  detailed  printing  occurs  again.     When 
IQ  reaches  three,    cycle  convergency  is  undoubtedly  satisfactory  and  line 
60  causes  the  program  to  proceed  to  the  next  element. 

Line  40002  is  the  end  of  the  cycle  or  JJ  iterative  DO  LOOP.      Line 
40000  ends  the  element  DO  LOOP.     Most  of  the  remainder  of  the  sub- 
routine concerns  printing.      There  are  however  a  few  calculations  remaining, 

In  lines   760  through  40001  the  work  done  in  all  elements  is  summed 
as  is  the  heat  transfer  and  the  friction.      There  are  now  two  types  of 
kinetic  energy  that  have  been  calculated.      In  each  element  the  infinitesimals 
have  been  summed  to  find  the  total  dewar  bound  kinetic  energy,    ZKE(I),    in 
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element  I.     The  kinetic  energies  in  either  the  N'th  or  M'th  element,    ZKE(N) 
or  ZKE(M),    could  serve  as  an  approximation  of  the  kinetic  energy  ejected 
into  the  dewar. 

The  kinetic  energy  has  also  been  calculated  for  each  time  interval 
in  line  550  and  summed  up  in  line  78  0  as  SUMEKE.      The  quantities  ZKE(N), 
ZKE(M),    and  SUMEKE  are  all  printed  and  found  to  agree  closely. 

The  purpose  of  the  program  is  to  compare  the  sum  of  work  done 
by  all  the  elements  with  the  kinetic  energy  ejected  into  the  dewar.      The 
best  comparison  is  a  ratio  of  the  kinetic  energy  to  the  work.      Lines  785 
and  787  compute  three  such  ratios.      The  quantity  EKEOSW  is  SUMEKE/ 
SUM  WORK;  ZKEMOSW  is  ZKE(M)/SUM  WORK  and  ZKENOSW  is  ZKE(N)/ 
SUM  WORK.     Because  it  is  based  on  the  most  reliable  method  the  quantity 
EKEOSW  is  used. 

9.     Detail  of  Subroutine  READ  DATA  in  Frequency  Calculations 

The  subroutine  READ  DATA  for  frequency  calculations  is  very 
similar  to  the  one  for  amplitude.      The  difference  is  that  much  of  the  input 
needed  for  amplitude  is  unnecessary  for  frequency.      The  input  quantities 
have  the  same  meaning  in  each  case.     Subroutine  READ  DATA  for  frequency 
is  given  in  Appendix  I. 

10.     Detail  of  Subroutine  ELEMENT  in  Frequency  Calculations 

The  subroutine  ELEMENT,    when  used  with  frequency  calculations, 
divides  the  duct  into  N  elements.      The  length  of  each  element  is  adjusted 
so  that  it  takes  sound  the  same  length  of  time  to  cross  each.      Thus,   the 
length  of  each  element  divided  by  its  velocity  of  sound  is  the  same,    and  the 
quantity,    DXXOVC,   is  the  same  for  each  element.      The  same  quantities 
are  computed  and  printed  for  each  element  as  they  were  for  amplitude 
calculations.     Additionally  the  velocity  of  sound  for  each  element,    C(I), 
is  computed  and  printed. 
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Line  10  sets  the  temperature  of  the  warm  end  boundary  of  element 
1,    TEMP(l),   to  the  input  value  of  the  w~rm  end  of  the  pipe,    TMP(l).     It 
also  sets  the  position  of  this  boundary,   XX(1),    equal  to  the  input  position 
of  the  warm  end  of  the  pipe,    S(l).      In  addition  line   10  computes  the  value 
of  NP1  =  N+l,    which  is  needed  later. 

Line  20  is  an  estimate  made  of  the  quantity  DXXOVC.     This  esti- 
mate is  the  length  of  the  whole  pipe,    S(32)  -  S(l),    divided  by  the  number 
of  elements,   N,    and  also  divided  by  the  velocity  of  sound  in  the  median 
element,   SQRTF(XK*R*TMP(16)). 

Lines  24  through  40002,   for  each  element  in  turn,    calculate  the 
sonic  velocity  in  an  iterative  routine,    and  determine  the  boundary  positions 
of  the  element.      Then  lines  47  through  49  test  to  see  if  the  cold  end  of  the 
N'th  element  is  in  its  proper  position.     If  not  an  adjustment  is  made  to 
the  quantity,    DXXOVC,    «-nd  the  program  returns  to  line  24. 


Thus,    line  25  estimates  the  sonic  velocity  of  the  element,    C(I), 
using  the  warm  side  temperature  TEMP(I).      Line  28  calculates  the  posi- 
tion of  the  cold  side.      Lines  29  through  33  calculate  the  temperature  of 
the  cold  side.     In  line  37,   the  average  temperature  of  the  element  is 
computed.      Then  the  tentative  sonic  velocity  of  the  element,    CCC,   is 
computed  using  the  average  element  temperature.      Line  38  tests  to  see 
if  the  tentative  and  estimated  velocities  are  nearly  the  same.     If  they  are 
not,    line  39  re -estimates  the  sonic  velocity  to  be  the  tentative  value  and 
the  program  returns  to  line  27.      If  they  are  nearly  the  same,    line  43 
establishes  the  sonic  velocity  as  being  the  tentative  value  and  the  program 
proceeds  to  the  next  element. 

At  line  47  the  positions  of  all  the  elements  have  been  computed. 
In  line  47  a  test  is  made  to  see  if  the  cold  side  position  of  the  N'th  element 
closely  coincides  with  the  position  of  the  dewar  end  of  the  pipe.     If  this  is 
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not  sufficiently  close,    the  quantity,    DXXOVC,    is  adjusted  in  line  49  and 
the  program  returns  to  line  24  repetitively  until  it  is  sufficiently  close 
after  a  number  of  iterations.      When  it  is  sufficiently  close  the  program 
goes  to  line  50  and  a  number  of  quantities  are  computed  and  printed  in 
a  manner  similar  to  that  of  subroutine  ELEMENT  in  amplitude  calculations. 
In  addition, three  other  quantities  are  calculated.     In  the  foregoing  iteration 
the  velocity  of  sound, C(I), was  determined  for  each  element.     In  line  55, 
the  parameter,  Z(I),  is  calculated.     This  is  simply  the  product  of  the  sonic 
velocity  and  the  density  of  the  gas  in  the  element,    RHO(I).     Also  computed 
is  its  reciprocal,    A(I).      These  two  parameters  are  needed  in  pressure  and 
velocity  calculations  which  determine  frequency.      The  remainder  of 
ELEMENT  is  for  print  out  purposes. 

11.     Detail  of  Program  FREQENCY 

The  program  computes  the  frequency  of  thermal  oscillations  from 
the  known  temperature  gradient  in  the  pipe,    using  equations  derived  in 
Appendix  A.     In  the  theory,    a  sudden  overpressure  is  assumed  in  the 
dewar  which  travels  towards  the  closed  end  according  to  the  velocity  of 
sound  and  the  derived  equations.     In  the  program  of  this  report  the  equiva- 
lent assumption  is  made  that  there  is  a  sudden  underpressure  in  the  pipe. 
Frequencies  calculated  by  the  two  methods  are  identical;    the  usage  is  one 
of  convenience.     In  either  case  the  initial  gas  velocities  are  zero. 

The  time  required  for  sound  to  traverse  one  element,    DXXOVC, 
is    one    computation  time.     At  the  start  of  the  program,    with  everything 
at  rest,    the  time,    J,    is  considered  to  be  one.     After  the  passage  of  one  time 
interval,    DXXOVC,    J  becomes  two.     At  the  time  [  J  =  2]    a  pulse  has  tra- 
versed from  the  dewar  across  element  N.      At  the  time  [  J  =  l]   the  program 
inserts  initial  values  of  overpressure,    P,    and  velocity,    U.     At  the  time 
[  J  =  2]   the  only  change  has  occurred  in  element  N,  and  the  program  deter- 
mines   the  proper  values.     After  these  first  two  preliminary  calculations, 
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the  program  computes,    at  each  time  J,    the  overpressure  and  velocity 
in  each  element.     Because  of  computer  storage  limitations,    the  maximum 
value  of  J  is  102.      Then,   if  desired,    the  pressures  and  velocities  in  the 
N  elements  are  printed  for  the  values  of  J  from  1  to  100.      The  values  of 
pressure  and  velocity  in  each  element  at  times  [  J  =  l]    and  [  J  =  2]    are 
now  set  equal  to  the  values  at  [J=10l]    and  [  J  =  102],    and  the  calculations 
proceed  for  another  hundred  time  intervals.     On  each  occasion  that  J  is 
thus  shifted,   the  quantity  KOUNT  is  increased  by  100  and  the  quantity  J2 
is  always  equal  to  J+KOUNT,   with  the  result  that  J2  is  the  number  of  the 
calculation  time.     Because  of  the  nature  of  the  equations  used,    a  different 
route  through  the  program  is  used  for  times  when  J  is  even  or  odd. 

When  first  devised,   the  frequency  was  determined  by  plotting  the 
closed  end  pressure  versus  time.      Then  the  number  of  cycles  per  second 
was  determined  by  examination  of  the  plot.     In  the  program  presented  here, 
the  frequency  is  determined  in  the  program  itself.     Because  of  the  some- 
what erratic  nature  of  the  overpressure,    the  routine  used  to  calculate  fre- 
quency is  quite  complicated. 

In  principle  the  method  is  as  follows.     During  each  calculation  time 
for  each  element  the  quantity,  PMIN( I,  JK),  is  determined.     If  the  overpres- 
sure,   P(I,  J),    at  the  time  is  less  than  the  previously  determined  PMIN(I,  JK), 
then  the  value  of  PMIN(I,  JK)     is  adjusted  downward  to  this  new  value. 
Simultaneously  the  quantity,    TMIN(I,  JK),    is  set  equal  to  the  calculation 
time,    J2,    at  which  the  overpressure  went  to  a  new  minimum.      The  sub- 
script,   I,    refers  to  the  element  number,    and  the  subscript,    JK,   to  the 
cycle  number  that  is  in  progress.     Initially  the  cycle  number  and  JK  are 
unity  and  the  closed  end  overpressure  is  -1000.   The  closed  end  pressure 
then  increases  to  a  maximum  overpressure  and  decreases.     Eventually 
it  becomes  negative.     At  this  time  a  new  low  is  being  approached  for  a 
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new  cycle  and  JK  is  increased  by  one.     When  a  new  minimum  pressure, 
PMIN(I,    JK),    occurs  at  a  new  time,    TMIN(I,    JK),    the  process  is  repeated. 

Without  some  safety  provision,   the  cycle  number,    JK,   would  in- 
crease prematurely.     This  results  from  the  fact  that  during  a  pressure 
trend  that  is   generally  up  or  down,    there  may  well  be  a  number  of  tem- 
porary reversals.     This  reversal  occurring  just  as  the  overpressure  cros- 
sed zero  would  erroneously  increase  JK.     The  safety  provision  introduced 
is  the  quantity,    MAXCOUNT.     At  each  change  of  JK,    MAXCOUNT  is  set 
to  zero.     During  each  calculation,    if  the  overpressure  is  slightly  over 
zero,    one  is  added  to  MAXCOUNT.      Until  MAXCOUNT  is  greater  than 
N,    the  quantity  JK  cannot  change.      By  this  time  the  danger  of  inadvertently 
having  the  overpressure  cross  zero  is  past. 

TMIN  of  the  first  element  is  used  as  a  criterion  of  the  interval  of 
one  cycle.     This  quantity  and  others  are  computed  for  all  the  elements. 
These  are, 


TMAX(I,  JK 
PMAX(I,  JK 
UMAX  (I,  JK 
TMIN(I,  JK) 
PMIN(I,  JK) 
UMIN(I,    JK) 

The  interval 


,   the  calculation  time  of  the  maximum  overpressure 
,   the  maximum  overpressure  of  cycle  JK 
,    the  maximum  velocity  in  cycle  JK 
the  calculation  time  of  the  minimum  overpressure 
the  minimum  overpressure  of  cycle  JK 
the  minimum  velocity  in  cycle  JK 


ZINTERVL,    of  one  cycle  is  obtained  by  subtracting 
a  given  value  of  TMIN  from  the  next  value  of  TMIN.     The  period,    PERIOD, 
is  this  quantity  multiplied  by  the  calculation  time,    DXXOVC.      The  recipro- 
cal of  any  one  period  is  the  frequency,    FREQ,    and  the  average  frequency, 
AVG  FREQ,    is  the  average  of  the  first  nine  frequencies  calculated.     The 
average  frequency  is  then  used  in  the  subroutine  AMPLITUD.      Program 
FREQENCY  is   given  in  Appendix  K. 
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Lines  5  and  10  call  the  subroutines  READ  DATA  and  ELEMENT. 

Line  20  sets  the  initial  value  of  KOUNT  at  -  100  so  that  when  the 
calculation  starts  and  100  is  added,    the  value  of  KOUNT  during  the  first 
100  calculations  will  be  zero.    The  quantity  SAMBDA  is  the  value  of  the 
initial  assumed  under   pressure  in  the  pipe.     GAMBDA  is  a  small  quantity 
used  as  a  criterion  to  determine  MAXCOUNT. 

Line  25  sets  the  values  of  several  quantities  which  are  used  later. 

Lines  30  and  40000  sets  the  pressure  in  all  elements  at  times  [  J  =  1] 
and  [  J  =  2]    to  -  SAMBDA  and  the  velocities  to  zero. 

Lines  40000  and  3  5  set  the  initial  values  of  quantities  which  may 
change  later. 

The  first  cycle  starts  at  its  minimum;  therefore,    line  40  sets 
PMIN  equal  to  -  SAMBDA.     For  other  cycles  it  is  set  to  +  SAMBDA  to 
insure  the  low  value  calculated  later  is  the  proper  one. 

Line  50  sets  the  velocity  and  pressure  in  the  N'th  element  and 
second  calculation  time  equal  to  initial  values. 

In  line  60  the  value  of  KOUNT  is  increased  to  zero  before  the  first 
series  of  a  hundred  calculations,    and  increased  by  100  after  each  series. 

Lines  65  through  141   are  the  main  DO  LOOP  that  calculates  the 
pressures  and  velocities.     In  line  65,    if  the  calculation  time  is  the  3rd, 
or  any  other  odd  time,    the  program  goes  to  line  101.      If  even,    it  goes  to 
line  121. 

Lines  101  through  115  are  the  DO  LOOP  that  calculates  odd  time 
overpressure,   velocity,    and  parameters  associated  with  PMAX  and  PMIN. 
The  overpressures  and  velocities  are  calculated  in  lines  105  through  113 
by  formulas  derived  in  Appendix  A.      Lines   2114  through  2129  determine  if 
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the  overpressures  and  velocities  calculated  are  a  maximum,    a  minimum 
or  neither  for  the  cycle  and  set  up  a  new  cycle  when  MAXCOUNT  is  greater 
than  N  and  the  closed  end  overpressure,    P(l,  J)  is  negative. 

Lines  121  through  5141  perform  the  corresponding  calculations  for 
even  values  of  J. 

At  line  150  the  overpressures  and  velocities  have  just  been  calcu- 
lated by  the  DO  LOOP  ending  in  line  141  for  a  hundred  values  of  J.     In  line 
150,   if  KOUNT  equals  zero  or  is  greater  than  1000,    detailed  printing  of 
pressures  and  velocities  is  desired  and  carried  out  by  lines  160  through 
225. 

In  line  226,   if  not  enough  cycles  have  been  completed,   the  program 
jumps  to  line  229  where  the  overpressures  and  velocities  of  the  101st  and 
102nd  calculations  times  are  transferred  to  the  1st  and  2nd  time,  and 
control  returns  to  line  60  for  computation  of  another  hundred  values  of  J. 

In  line  226,   if  enough  cycles  have  been  completed,    the  program 
continues  to  line  6227  for  final  calculations. 

Line  6228  finds  the  number  of  calculation  times,    ZINTERVL(J),   for 
each  cycle. 

Line  6229  converts  ZINTERVL(J)  to  cycle  period,     PERIOD(J),   by 
multiplying  by  the  calculation  time,   DXXOVC. 

Line  6230  takes  the  reciprocal  of  PERIOD(J)  to  find  the  frequency 
of  each  cycle,   FREQ(J). 

Line  6231   adds  all  nine  frequencies  together. 

Line  6232  determines  the  final  answer  needed  for  AMPLITUD,   which 
is  the  average  frequency,   AVGFREQ. 

The  remainder  of  the  program  is  for  print-out. 
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12.     Results  and  Conclusions 

Typical  results  from  the  frequency  program  are  shown  in  figure 
3.      To  obtain  this  plot  a  version  of  the    program  is  used  which  does  not 
compute  the  frequency.     Instead,    the  closed  end  overpressures  are  print- 
ed and  then  plotted.     The  plot  shows  the  almost  sinusoidal  overpressure 
obtained.     It  also  shows  the  tendency  for  the  calculated  overpressure  to 
be  erratic  even  though  the  general  trend  is  clearly  recognizable.     Since 
the  temperature  profile  for  the  calculations  is  from  figure  4  which  is  from 
the  work  of  Bannister  (1965),    it  is  possible  to  compare  measured  with 
calculated  frequencies. 

The  curve  of  figure  3  shows  three  oscillation  cycles  over  a  range 
of  321  calculation  points  which  followed  3984  previous  calculation  points. 
Over  the  whole  calculation,    the  period  averaged  107.2  calculation  inter- 
vals of  0.  000302  seconds  each.     From  this  the  calculated  frequency  is  30.9 
cps.     From  figure  3,   the  period  is  about  105  calculation  intervals,    or 
0.  0317  seconds,    and  the  frequency  is  31.5  cps.     From  the  program  of 
this  report  which  averaged  the  first  10  cycles,    the  frequency  is  31.  25  cps. 

It  should  be  noted  that  figure  3  is  from  a  slightly  different  program 

than  the  one  described  in  detail  in  this  report.     In  figure  3  an  initial  over- 

2 
pressure  in  the  dewar  of  10130  dynes/cm     is  assumed.     In  the  described 

2 
program  an  initial  underpressure  of  1000  dynes/cm     is  assumed  through- 
out the  pipe.     This  difference  has  no  effect  on  the  calculated  frequency. 

The  calculated  frequency  is  needed  in  the  determination  of  pres- 
sure amplitude.     In  all  the  amplitude  determinations  of  this  report  the 
temperature  profile  of  figure  4  is  used,    so  the  frequency  of  each  calcu- 
lated point  is  taken  as   31.  25  cps.      The  program  is   designed  however  to 
accommodate  any  predescribed  temperature  profile  by  modifying  input  data. 
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The  amplitude  calculations  depend  on  heat  transfer  coefficients 
and  friction  factors.     Neither  of  these  quantities  is  precisely  known  for 
the  type  of  fluid  movement  that  exists  in  thermal  oscillations.     In  Mc Adams 
(1954)  is  information  for  fully  developed  flow.      Unfortunately  the  flow  in 
ftiis  case  is  far  from  fully  developed.     Furthermore,    it  is  intended  that  a 
high  value  of  friction  factor  is  to  be  used,    to  allow  for  inlet  turbulence. 

Accordingly,    the  results  are  presented  for  two  types  of  heat  trans- 
fer coefficients.     With  the  one  type,    heat  transfer  coefficients  are  obtain- 
ed directly  from  the  text  book.     With  the  other,    they  are  multiplied  by  0.9. 

Two  methods  are  used  for  presenting  friction  factor.     With  one 
method,    friction  factors  are  assumed  to  be  constant  during  a  run  regard- 
less of  the  momentary  Reynold*s  number.      Then  a  series  of  runs  are 
presented,    each  with  a  given  constant  friction  factor.     With  the  other 
method,    friction  factors  are  calculated  from  a  text  book  formula  for  each 
momentary  Reynold's  number,    but  the  calculated  value  is  increased  by  a 
constant  multiplier  for  each  run.     Then  a  series  of  runs  is  presented, 
each  with  a  constant  friction  factor  multiplier. 

Three  figures  are  presented,   using  a  variety  of  friction  factors 
and  two  coefficient  multipliers,    showing  the  effect  of  pipe  diameter  on 
the  amplitude  of  oscillations.     Figure   5  uses  unmodified  text  book  heat 
transfer  coefficients  and  a  series  of  four  constant  friction  factors.     Fig- 
ure 6  uses  unmodified  text  book  heat  transfer  coefficients      and  a  series 
of  five  friction  factor  multipliers.     Figure  7  uses  text  book  heat  transfer 
coefficients  times  a  multiplier  of  0.9     and  the  same  series  of  five  friction 
factor  multipliers. 

The  general  shape  of  all  the  curves  is  the  same.     At  the  small 
diameter  end,   the  oscillations  tend  towards  zero  amplitude.     As  the 
diameter    increases,   the  calculated  amplitude  increases  to  a  maximum, 
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and  then  decrease  asymptotically  towards  zero  at  large  diameters.     As 
might  be  expected,    the  higher  friction  factors  result  in  lower  peak    over- 
pressures in  figures   5,    6,    and  7.     In  figures  6  and  7  the  effect  of  heat 
transfer  coefficients  may  be  seen.     In  all  cases,    reducing  the  heat  trans- 
fer coefficient  multiplier  from  1„  0  to  0.9  reduced  the  peak  overpressure. 
As  an  approximation,    the  reduction  is  about  30  percent. 

It  is  of  interest  to  compare  the  calculatjd  overpressures  with  the 

data  from  Bannister  (1965).     For  a  0.  63  5  cm  (1/4  inch)  diameter  pipe 

2 
the  measured  peak  to  peak  pressure  oscillation  was  220,  000  dynes/cm 

2 
(3.  19  psi)  which  corresponds  to  a  peak  overpressure  of  110,  000  dynes/cm    . 

2 
On  figure  5  the  constant  friction  factor  that  produces   110,  000  dynes/cm 

is  a  value  of  0.  0065.     On  figure  6,    by  interpolation,    the  friction  factor 
multiplier,   which  produces  this  overpressure,  is  about  1.  57,    and  on  fig- 
ure 7,    by  extrapolation,    tne  value  is   1.37. 

The  constant  friction  factor  method  is  not  favored.     At  small  dia- 
meters,   the  Reynold's  number  will  be  smaller  and  the  friction  factor 
larger  than  at  large  diameters.     As  a  result,   the  maximum  overpressure 
occurs  at  too  low  a  diameter  as  is  the  case  on  figure  5. 

It  is  not  usually  expected  that  heat  transfer  coefficients  are  known 
precisely.     Therefore  the  effect  of  diameter  on  amplitude  can  be  followed 
equally  well  with  a  heat  transfer  coefficient  multiplier  of  1.0  or  0.  9.     The 
value  of  1.  0  from  figure  6  is  therefore  used.     From  the  one  datum  point 
the  friction  factor  multiplier  of  1.  57  has  been  selected.     From  the  figure 
it  is  found  that  the  maximum  overpressure,    for  the  given  temperature 
gradient,    occurs  at  0.  3175  cm  (1/8  inch)  and  its  value  is  about  480,  000 
dynes/cm     (7  psi),    and  double  that  value  for  peak  to  peak  variation.     At 
0.  63  5  cm  (1/4  inch)  the  value  coincides  with  test  data  and  at  large  dia- 
meters tends  towards  zero. 
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Considering  that  no  allowance  exists  in  the  program  for  entrance 
turbulence,    and  further  considering  that  friction  factor  determination 
cannot  be  precise,    it  is  quite  reasonable  to  have  to  multiply  by  1.  37  or 
1.  57  to  correlate  calculated  oscillation  overpressures  with  test  data. 

It  is  concluded  that  the  program  in  its  present  form  calculates 
overpressure  peak  amplitudes  that  are  approximately  correct.     Therefore, 
the  qualitative  explanation  of  the  phenomenon,    given  in  this  report,    which 
generally  follows  Rayleigh  (1945),    is  correct.     Further, the  heating  effect 
which  evaporates  helium  is  the  result  of  the  operation  of  the  work  cycles 
which  expel  kinetic  energy  and  warmed  gas  into  the  liquid. 

Thus  the  method  of  oscillation  suppression,  by  inserting  a  knotted 
string  in  the  pipe     is  in  effect  increasing  the  coefficient  of  friction. 
Similarly  the  insertion  of  glass  wool  or  other    packing  does  the  same  thing. 
Clement  and  Gaffney  (1954)  found  that  a  vacuum  insulated  line  does  not 
oscillate.     This  case  corresponds  to  a  heat  transfer  coefficient  multiplier 
of  about  zero.     Since  a     0.  9  multiplier  greatly  reduces  oscillations,    it  is 
expected  that  a  multiplier  close  to  the  value  of  zero  would  essentially 
eliminate  them.     Furthermore,    since  the  oscillations  are  a  thermal  cycle, 
a  source  of  heat  is  needed  to  provide  the  cycle  work.     Complete  insulation 
eliminates  the  source.     In  the  same  reference  it  is  found  that  closure  of 
the  cold  end  suppresses  oscillations.     For  this  case  there  is  no    kinetic 
energy  ejected  into  the  dewar.     Any  chance  oscillations  that  start  come  to 
equilibrium  when  the  sum  of  all  the  work  terms  equal  the  kinetic  energy. 
This  equilibrium  occurs  when  the  work  terms  are  all  zero,    or  the  oscil- 
lations are  non-existent. 

Although  the  program  does  simulate  frequency  and  amplitude  of 
thermal  oscillations,    it  has  a  defect  which  should  be  improved  upon.     The 
assumption  that  the  maximum  overpressure  varies  as  a  cosine  function 
from  one  end  to  the  other  of  the  pipe  is  not  needed.     The  assumption  that 
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the  closed  end  pressure  is  sinusoidal  with  time  is  correct  and  should  be 
retained.     Then  using  the  gas  laws  and  Newton's  laws  of  motion  the 
boundaries  of  the  elements  can  be  calculated  one  by  one  from  the  closed 
to  the  open  end.     Simultaneously  the  instantaneous  pressure  variation  at 
all  of  the  boundaries  can  be  obtained.     Such  a  calculation  would  be  an 
improvement  over  the  present  one  and  is  recommended  for  future  con- 
siderations.    Additional  work  is  needed,    particularly  in  experimental 
phases  such  that  more  appropriate  comparisons  can  be  made  with  the 
mathematical  model. 

Although  calculations  have  been  performed  and  comparisons 
made  for  one  particular  case,   the  program  is  expressly  designed  to 
accommodate  different  fluids,    tube  lengths  and  temperature  gradients 
along  the  tube. 
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14.     Appendix  A 

To  analyze  the  behavior  of  pressure  pulses  traversing  a  non 
homogeneous  medium  it  was  necessary  to  make  an  assumption.      The  gas 
temperature  in  the  tube  varies  continuously  from  one  end  to  the  other. 
If  an  attempt  is  made  to  analyze  this  situation,    an  infinite  number  of 
reflected  pressure  pulses  would  be  obtained,    which  cannot  be  handled  in 
a  satisfactory  manner.     It  therefore  was  assumed  that  the  behavior  of  the 
continuous  temperature  gradient  could  be  closely  approximated  by  an 
equivalent  gradient  composed  of  a  large  number  of  constant  temperature 
elements.     As  the  number  of  these  elements  becomes  sufficiently  large, 
the  behavior  of  the  continuous  gradient  and  the  discontinuous  gradient 
should  become  identical.     Such  an  assumption  is  quite  common  in 
analyses  of  physical  situations.      The  assumption  makes  no  requirement 
of  the  length  of  each  gas  element  other  than  the  tacit  one  that  each  element 
length  should  be  small.     As  a  convenience  in  analysis,    the  element  lengths 
were  selected  such  that  the  time  required  for  sound  to  traverse  each  gas 
element  would  be  the  same.     As  a  further  convenience  in  analysis  an  even 
number  of  elements  was  used. 

In  the  derivations,    all  pressures  were  considered  as  absolute. 
However,    since  the  derived  expressions  showed  pressure  differences,    it 
was  often  possible  in  the  digital  computer  program  to  ignore  the  absolute 
pressure  level  and  to  consider  only  over  or  underpressure.     Of  course, 
in  determining  gas  densities,    the  absolute  pressure  must  be  used.      The 
quantity,    J,    is  used  in  analysis.     When  J  is  one  the  overpressure  has  just 
occurred.     After  the  length  of  time  it  takes  sound  to  traverse  one  element, 
J  has  increased  to  two.     After  each  such  time  increment, J     is  increased  by 
one. 

At  time  [J=  l]  ,   with  reference  to  figure  8,    the    process  is  just 
ready  to  start.     In  each  element  the  overpressure  and  gas  velocity  are 
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zero.     In  the  dewar  there  is  a  small  overpressure  X.     A  moving  sonic 
front  exists  at  the  tube  end  between  the  dewar  overpressure  and  the  tube. 

At  time  [  J  =  2]   the  dewar  overpressure  has  penetrated  across  the 
gas  element  nearest  the  dewar,    element  [I  =N]  .      Because  of  the  dis- 
continuity in  gas  properties  between  element  [  I  =N]    and  [  I  =  N  -  l]  ,    a 
portion  of  the  overpressure  is  transmitted  into  element  [I  =  N  -  l]    and 
a  portion  is  reflected  through  element  [  I  =  N]  .     The  values  of  the  result- 
ing gas  velocities  and  pressures  will  be  derived  later.     At  time  [  J  =  2] 
there  is  no  longer  a  sonic  front  at  the  tube  end  but  there  now  is  one  be- 
tween elements  [  I  =  N]    and  [  I  =  N  -  l]  . 

At  time  [  J  =  3]   the  transmitted  wave  has  traversed  element 
[  I  =  N  -  1]    and  is  at  the  border  of  elements  [  I  =  N  -  l]    and  [  I  =  N  -  2]  . 
At  the  same  time  the  reflected  wave  is  at  the  tube  end  between  the  dewar 
and  element  [  I  =  N]  .     At  time  [  J  =  3]   there  is  a  sonic  front  between  the 
dewar  and  element  [  I  =  N]  ,    and  between  elements  [  I  =  N  -  l]    and  [  N  -2]  . 

At  each  succeeding  time  there  will  be  a  new  pattern  of  sonic  fronts. 
After  the  transmitted  pulses  have  reached  the  closed  end  it  may  be  seen 
that  the  pattern  will  be  continuously  repetitive.     At  each  time  [  J  =  odd  num- 
ber]  there  will  be  sonic  fronts  at  the  two  tube  ends  and  at  the  dewar  end 
of  all  even  gas  elements.     (This  was  the  reason  for  using  an  even  number 
of  elements.  )    At  each  time  [  J  =  even  number]  ,    there  will  be  no  sonic 
fronts  at  the  tube  ends,    but  there  will  be  some  at  the  dewar  end  of  all  odd 
gas  elements. 

For  a  tube  with  N  =  8  gas  elements  this  situation  is  illustrated  on 
figure  8  . 

From  figure  8   it  may  be  seen  that  at  any  time  [  J  =  odd]    all  gas 
elements  are  grouped  in  pairs,    with  each  pair  surrounded  by  sonic  fronts. 
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For  any  time  [  J  =  even]    most  of  the  gas  elements  are  grouped  in  such 
pairs,    but  elements,    [  I  =  l]    and  [  I  =  N]  ,    are  arranged  independently. 

In  the  following  derivations  the  gas  velocities  and  pressures  are 
derived  in  terms  of  gas  velocities  and  pressures  at  previous  times.     The 
first  derivation  is  for  element  pairs  surrounded  by  sonic  fronts.     Separate 
derivations  then  follow  for  the  cases,    [  Time  J  =  even,    Element  1=1  and 

I  =  N]. 

On  figure  9  at  time  [  J  -2]   there  are  two  regions,  [  I]    and  [I  -  1], 
surrounded  by  sonic  fronts.      Later  at  time  [  J]   they  are  again  in  the  same 
condition.     At  time  [  J  -  2]   two  pressure  pulses  have  just  started  towards 
the  center.     At  time  [  J  -  1]   they  have  met  at  a  density  discontinuity,    and 
at  time  [  J]   they  are  widely  separated  again.     At  the  intermediate  times 
[t   ]    and  [t   ]   the  pressures,    velocities  and  other  gas  functions  9  are 
known  on  each  side  of  the  sonic  fronts  in  terms  of  functions  cp  at  times 
[J  -2,    J  -  1,    and  J]  . 

At  time  [t   ]    on  the  left  side,    from  continuity 

dM=(CI-l.    t-UI-l,    J-2)PI-1,    J-2   Adt  (1) 


dM=(CI-l,.1-U.-I,    J-1)PI-1,    J-l    Adt-  (2> 

From  momentum  it  follows  that 


(pi-i,  j.r  pi-i,  j-2>  = dM  <ui.i,  j.r  ui-i,  j-2>/<Adt)  •    <3' 

Combining  (1)  and  (3)  results  in 


(P  -   P  )  = 

V    1-1,    J-l         1-1,    J-2; 

(CI-1,   tx~   UI-1,    J-2>  PI-1,    J-2  (UI-1,    J-l"  UI-1,    J-2>' 
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Figure  9.     Element  Pair  Surrounded  by  Sonic  Fronts. 
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Combining  (1)  and  (2)  leads  to 


(CI-l.t   -  UI-l,J-2>  PI-l,J-2  "  (CI-l>tl-  UI-1,  J-l>  "i-l.  J-l  •     <5> 


By  algebraic  manipulation  of  (5) 


(CI-1,   f  UI-1,    J-2>  = 


(6) 


(ui-i,  j-r  ui-i,  j-z)  pi-i,  j-i/(pi-i,  j-r  pi-i,  j-2>  • 

Solving  (4)  for  (C  -  U  )  and  multiplying  by  (6) 

1~  1  j     £  ..  x~  1 1     J  —  £* 

produces 

.2  =    PI-1,   J-l  '  PI-1,    J-2    PI-1,    J-l 
^    4l"      I"1'    J"2  "I-l.    J-l   "   "l-l.    J-2     'i-l,    J-2 


At  time  t     on  the  right  side,    by  the  same  operations 

(PI,J-2  "  PI,    J-l^  =  <CI,   tl+UI,    J-2>  PI,    J-2  <UI,    J-l"  UI,    J-2> 

(4)' 

(P  -  P  )    P 

^1         *'    J"2       "(PI,    J-2"PI,    J-l)     PI,    J-2  U 


At  time  t     on  the  left  side 
2 


(PI-1,    J-l"  PI-1,    j)  = 


(CI-1,   t2+UI-l,    J-1>PI-1,    J-1<UI-1,    J"UI-1,    J-l) 


(4)» 
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(C  +  U  T1)2=Pl-1'    J-'  "  "■•''   J   ^i        .  (7)" 

H.  t2      l-l,  J-i        p        j_x  -p        j   f         ^ 


At  time  t     on  the  right  side 


<pi,  j"  Pi,  j-i>  ■  <ci,  «2-  ui,  j-i>  \  j-i  <ui,  j"  ui,  j-i»      <4>'" 

P           -  P  p 

(C  -U  \2  =  _J±_i 1>   J"1         *>   J  (7)ih 

^z       *•   J"x  PI,   J"\  J-l     PI,  J-l       '  ' 

It  may  be  noted  that  where  the  velocity  of  the  fluid  is  small  com- 
pared with  sonic  velocity  and  where  the  density  in  a  given  element  changes 
only  slightly  with  the  passage  of  a  pulse,   the  velocity  of  sound  relative  to 
the  wall  is  the  usual  expression  for  sonic  velocity. 

Since  each  sonic  front  during  any  time  period  is  in  a  medium  of 
constant  properties,   each  front  is  at  a  constant  velocity.     From  this  fact 
and  from  continuity  it  follows  that  pressure  and  gas  velocity  between  two 
adjacent  sonic  fronts  is  the  same  on  each  side  of  the  density  discontinuity. 

PI-1,    J-2  =  PI,    J-2  *8) 

°I-1.   J-2  "  UI,   J-2  .  (9) 

P  =     P 

1-1,    J  I,    J  (10) 

ui-i,  j-uW  <"> 

For  small  pressure  pulses  the  equations  may  be  simplified.  The 
density  in  a  given  region  will  not  change  appreciably,  and  the  gas  veloc- 
ities will  be  small  compared  with  sonic  velocities.     After  this 
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simplification  and  using  (9)  through  (12),    (4)  becomes 

P  -P  =Cp(U  -U  )  (12) 

1-1,   J-l         I,    J-2         1-1  MI-1  |UI-1,   J-l        I,    J-2'  lx   ' 

PI,   J-2  "  PI,   J-l  "  CI  PI  (UI,    J-l  -  \   J-2>  <12>' 

PI-1,   J-l  "  PI,   J  "  C!-l  PI-1  (UI,   J-UI-1,   J-l'  <12)" 

PI,   J  "  PI,   J-l  "  CI  h  (UI,   J  "  UI,   J-l»    •  (12' 

Adding  (12)  to  (12)'  produces  the  intermediate  equation 

P  -P  =C  p(U  -U  J 

1-1,    J-l  I,    J-l  1-1,    HI-1  l    I-l.J-1        I,    J-2' 

+  C     p    (U  -  U  )  . 

T     I  MI  v    I,   J-l  I,    J-2; 

Adding  (12)"  to  (12)'"  gives 

P  -P  =Cp(U-U  ) 

1-1,    J-l        I,    J-l  1-1   HI-1  V    I,    J  1-1,   J-l1 

+  ci  pi  (ui,  j  "  ui,  j-i'  • 

Upon  adding  the  intermediate  equations 

<ui,  j  -  ui,  j-2'  ■ 2  \\  iXXl'1  ■  <13) 

Upon   subtracting  the  intermediate  equations 


1  \  /      1      \ 


ui.  J-i  ^c 

<ui,  j  +  ui,  j-2>  =  2 1  iV7    i      x — <13>' 

V/   ic.-i '.V 
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Similarly 

(pi,  j  -  pi,  j-1  = z  rt  vtt.  \  J'\  •  (14» 


PI  J-l  (CI-1  "l-l* +  PI-1  J-l  (CI  "i* 
,pw  +  pt  j-z>  ■  2  v^CmIi  -(14'' 

Equations  (13)  and  (14)  determine  the  gas  velocity  and  pressure 
at  any  time  [  J  =  odd]    for  any  gas  element  [i  =  even]    in  terms  of  previous- 
ly determined  velocities  and  pressures.     In  each  element  [  (I  -  1)  =  odd] , 
the  velocities  and  pressures  are  identical  to  those  in  the  elements  [I  = 
even]   because  of  (8),    (9),    (10),    and  (11). 

Equations  (13)  and  (14)  also  determine  the  gas  velocity  and  pres- 
sure at  any  time  [  J  =  even]   for  all  gas  elements  [I  =  odd]    except  for  the 
first  element  [  I  =  1]  .     In  each  element  [  (I  -  1)  =  even]  ,  the  velocities 
and  pressures  are  identical  to  those  in  the  elements  [I  =  odd]  . 

For  all  times  [  J  =  even]  the  pressure  in  the  element  [  I  =  N]   is 
dewar  pressure.     The  gas  velocities  are  obtained  with  the  aid  of  figure  10. 
Using  the  principles  of  momentum  and  continuity  as  before 

(P  -  P         ) 

TT  TT  V     N,    J-l  N,  J; 

N.   J"      N,   J-I-P^^^C^^+U^    j^)      ' 

For  small  pulses  the  density  does  not  change  appreciably  from  one  time 
to  another.     Also  the  velocity  of  the  gas  is  small  compared  with  sonic 
velocity.     So  for  small  pulses  for  all  times  [  J  =  even]   the  gas  velocity  in 
the  element  [  I  =  N]    may  be  found  from, 


61 


Element  N 


0> 


WjJ--/ 


dcwar 


Time   t/V  *  odcJ 


9*1,7-1 


N,J 


°N,t, 


T/rne     t, 


9m,j-/  I     9/i/r 


a 


%t 


T/'me     t,  feft 


0*, 


j 


T/'me    J~=  eifc/7 


Figure  10.     Sonic  Front  in  Element  N. 


62 


U  -U  _  PN,    J-l   '  PN,    J 

UN,    J       UN,    J-l"  PNCN  •  <l6> 

For  all  times  [  J  =  even]  ,    from  continuity,   the  velocity  in  the 
element  [I  =  l]    is  zero.     The  gas  pressures  are  obtained  with  the  aid  of 
figure  11.     Using  the  momentum  and  continuity  principles  again  and 
assuming  small  pressure  pulses  yields 

Pl,   J  "  Pl,   J-2  "  "  2  Cl  Pl  Ul,   J-l  =  "  2  Cl  Pl  U2,  J-l    •  <17> 

Equations  have  now  been  derived  which  show  gas  velocity  and 
pressure  in  each  of  N  gas  elements  at  any  calculation  time  J  provided  the 
velocities  and  pressures  are  known  for  the  previous  two  calculation  times, 
For  the  first  calculation  time  [  J  =  1] ,    all  velocities  and  pressures  are 
zero.     For  the  second  calculation  time  [  J  =  2]    all  velocities  and  pres- 
sures are  zero  except  for  the  element  [  I  =  N]  .     In  that  element  the  pres- 
sure is  ^  above  original  and  the  velocity  is  obtained  from  equation  (16). 


Nomenclature  for  Appendix  A 


Symbols 


2 
P  Pressure  dynes/cm 

C  Sonic  velocity  relative  to  tube  cm/sec 

3 
p  Density  g/crn 

U  Gas  velocity  (positive  toward  dewar)  cm/sec 

2 
X  Initial  overpressure  dynes /cm 

M  Mas  s  g 

cp  Sonic  function,   pressure,    density, 

velocity  as  determined  by  context. 

2 

A  Area  cm 
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Element  / 


0, 


,J~-2 


T?'me>      J~2  -  euzn 


Time      t/ 


0*>"     Time    tt+dt 


0. 


/.j--/ 


h 


*'     T/me  J-/  =catt 


'2.  J-/ 


T/'me     t> 


2<J~f   Time     tz+dt 


<Pi. 


j 


T/'me     J~=  ever? 


Figure  11.   Sonic  Front  in  Element  1, 
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Nomenclature  for  Appendix  A  (continued) 

Subscripts 

I  Gas  element,    I  =  1  at  closed  end,    N  at  open  end 

N  Gas  element  closest  to  dewar 

J  Calculation  time.     At  J  =  1  process  is  just  ready  to  start. 

Time  between  J  and  J  4*  1  is  time  required  for  sound  to 
traverse  one  gas  element. 
Where  double  subscripts  are  used,   the  first  refers  to  a 
gas  element,   the  second  to  the  calculation  time  interval. 
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1  5.    Appendix  B 

The  infinitesimal     net  work  done  by  a  gas  element  in  a  short  time 
interval  is, 

dW  =  P  dV  -  dF 


where 

dW 

= 

work  done             -   - 

ergs 

P 

= 

pressure                -   - 

dynes/cm 

dV 

= 

volume  change  -   - 

3 

cm 

dF 

= 

friction  work      -   - 

ergs     • 

From  the  gas  laws  we  have, 
PV       =    MRT 


dV 
V 

dT       dP 
T    "     P 

PdV 

=    MR  dT-VdP 

where 

R    = 

gas  constant 

T    = 

gas  temperature 

M  = 

gas  mass 

ergs/g-°K 
°K 
g  • 
The  temperature  rise  may  be  expressed  as, 

dQ  -  dW 


dT  = 


M  C 
v 


where 

dQ     =    heat  received         -   -   -      ergs 

C        =     specific  heat  -   -   -      ergs/g-°K. 

Substituting, 


PdV  =  ^?-   (dQ  -  dW)   -  VdP 

v 
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PdV  =  -^-  dQ  -  -^-  (PdV  -dF)  -  VdP 
v       v 


R        R 
—  dQ  +  —  dF  -  VdP 

v       v 
PdV  =  — 


V 


(v_i)      (Y-1)      VdP 
PdV  =  ^— J   dQ  +  ^—^   dF  -  -^=- 

y  y  y 


dW  =  PdV  -  dF  =  i^-^-J   dQ  - 
where, 


/VIl)dQ_VdP_dF 
Y    /  Y  Y 


Y  =    specific  heat  ratio. 

In  Fortran  notation  this  becomes, 
DWK      =    C9*DQ-VOLM/XK*(P(2)-P(l))-DFR/XK 
where, 

DWK      =      net  work  -  -  -       ergs 

XK  =      specific  heat  ratio 

C9  =      (XK-1)/XK 

DQ  =      heat  received  -  -  -      ergs 

3 
VOLM  =      average  volume  -  -   -       cm 

2 
P(2)       =      overpressure  at  end  of  interval  -   -  -       dynes/cm 

2 
P  (1)      =      overpressure  at  start  of  interval        -  -  -       dynes/cm 

DFR       =      friction  work ergs   . 
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16.    Appendix  C 

When  a  fluid  is  moving  in  a  pipe  there  will  be  a  frictional  force 
between  any  gas  element  and  the  wall.     This  is  expressed  in  equation 
(6-2c)  of  McAdams  (1954)  as 

T  =  £  v— y 

where 

t    =     force  per  unit  wall  area     -  -   -     dynes /cm 

f     =     Fanning  friction  factor 

3 
p     =      fluid  density  -  -  -     g/cm 

V   =      fluid  velocity  -   -   -     cm/sec 

From  this,   the  frictional  work  overcome  by  an  element  with  a 

given  area  is, 

-  PV   \ 
FdX^f  ii~  )  A  V  dt 


where 


where 


F       =     wall  force  -  -   -      dynes 

dX    =     distance  moved  -  -  -      cm 

2 
A      =     wall  contact  area       -  -  -      cm 

V       =     fluid  velocity  -   -   -      cm/sec 

dt      =     time  increment  -  -  -      sec 

In  Fortran  notation 

DFR  =  ABSF(FRF/2.  *AREA*ROE*VEL**3*DT) 


DFR      =     friction  work  -  -  -       ergs 

2 

AREA  =     element  wall  contact  area        -   -   -       cm 

3 

ROE       =     average  density  -   -  -        g/cm 

VEL       =     fluid  velocity  -  -   -       cm/sec 

DT  =     time  increment  -  -   -       sec 
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The  absolute  value  is  used  because  friction  work  is   always  a  positive 
number  to  be  subtracted  from  the  work  done. 
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17.     Appendix  D.     Program  THM  LOS 


PROGRAM  THMLOS 
CALL  READ  DATA 
CALL  ELEMENT 
CALL  AMPLITUD 
END  THMLOS 
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18.     Appendix  E.     Subroutine  READ  DATA  for  Amplitude 


SUBROUTINE  READ  DATA 
C  FRICTION  FACTOR  (FRF)  IS  READ»  BUT  IT  MAY  BE  CHANGED  IN  (HEAT  TR) 
C  IF  THE  INPUT  FRICTION  FACTOR  (FRF)  IS  ZER0»THE  QUANTITY  (IFRF)  IS  SET  EQUAL  TO 
C  ZERO  IN  SUBROUTINE  (AMPLITUO) ,  AND  THE  SUBROUTINE  (HEAT  TR)  COMPUTES  THE  HEAT 
C  TRANSFER  COEFFICIENT  ACCORDING  TO  A  FORMULA  IN  MC  ADAMS»  MODIFIED  By  THE 
C  MULTIPLIER  (NFR) . 
C 

C  IF  THE  INPUT  FRICTION  FACTOR (FRF) IS  NON  ZERO,  THE  QUANTITY  (IFRF)  IS  SET  EQUAL 
C  TO  UNITY  IN  SUBROUTINE  (AMPLITUD)  AND  THE  INPUT  VALUE  OF  THE  FRICTION  FACTOR 
C  IS  USED  REGAROLESS  OF  REYNOLDS  NUMBER  (REYNO) . 

COMMON/l/pSTART»DI AM, FRF, FACTOR, AVGFREQ,PM,N,NJJ,PI,AF,XK,CPG»M, 
1  TMP(32)fS(32)»TMRVIS(l6) »VIS (16) ,TCONGAS ( 12) » CONGAS ( 12) ,ZNU<16), 
1  GRPR<16) »CVG,RHODEW,TDEW,R,jD 
C0MM0N/4/NFR 
TYPE  REAL  NFR 
5   READ800ltPSTART,DI AM, FRF, FACTOR, AVGFREQ,PM,N,NJJ. NFR, JD 
8001   FORMAT(6F10.0»2I5»F5.0»I5) 

10   READ  8000f  (TMP(I) ♦1*1,32)  $  READ  8000 » <S < I ) » 1=1 ,32) 
15   READ  8000,  (TMR  VIS ( IN) ♦ IN=1 ,  16)  SREAD  8000, (VIS     (IN) , IN=1 , 16) 
20  READ8000HT  CON  GAS  ( IN)  ,  IN«1 ,  12)  $READ8000»  (CON  GAS  ( IN)  ,  IN*1 ,  12) 
25   rEAD  8000* (Z  NU    ( IN) , IN»1 , 16)  SREAD  8000,  (Gr  Pr   ( IN) , IN  =  1 . 16) 
8000   FORMATU6F5.0) 
CONVERSION  OF  INPUT  DATA  TO  UNITS  USED  IN  CALCULATIONS 

30   DO  40012  K*1.32  $S (K) =S (K) *2.54 
40012   CONTINUE   $DIAM=DIAM*2.54 

35    DO  40011  Kxl,16  $vIS <K) *yIS <K) /I  000  000. 
40011   CONTINUE 

40    DO  40010  K=l,12  SCONGAS  (K) .CONGAS (K) •10  000. 
40010   CONTINUE 

C  VISCOCITY  OF  GAS»VlS»IS  REAO  IN  MICROPORE  AND  CONVERTED  TO  POISE 
C  C0NDUCTIVITY»C0N  GAS,  IS  READ  IN  MILLIWATTS/ (CM  DEG  K)AND  CONVERTED  TO 
C  ERG/(S  CM  DEG  K) ,  TMR  VIS, DEG  K.GOES  WITH  VIS 
C  T  CON  GAS,  DEG  K,  GOES  WITH  CON  GAS,    FOR  THE  GRASHOF 

C  PRANDTL  PRODUCT.GR  PR,VERSUS  THE  NUSSELT  NUMBER, Z  NU,  THE  LOG  10  OF  THE 
C  QUANTITIES  ARE  READ  IN  WITH  NO  CONVERSION.   THE  LENGTH  S  IS  READ  IN 
C  INCHES  AND  CONVERTED  TO  CENTIMETERS.  THE  SAME  IS  TRUE  FOR  DIAM 
C 

45   WT=4,003SRU=8.317E07$R3RU/WT$PIs3.l4l5926536$AF=PI/4.«DIAM*«2 
50   XK=5./3.$  CPG=R*XK/(XK-1.)$CVG«R/(XK-1.) 
55   T  DEWs4.3  $M=N*N/2  SRHODEw»PSTART/R/TDEW 
COMPUTATIONAL  VALUES  ARE  PRINTED  EXCEPT  VISCOCITY  IS  PRINTED  IN  MICROPOISE 
60   PRINT  9310 
9310   FORMAT (//) 

65   PRINT  9000.PSTART, DIAM, FACTOR, AVGFREQ,PM,N,NJJ, NFR, JD 

70   PRINT9001,(TMR  VIS (K) ,K-1 , 16)  SPRINT  9003» (VIS  (K) ,K»1 , 16) 

75   PRINT  9002»(T  CON  GAS (K) ,K»1 ,12)  SPRINT  9004, (CON  GAS (K) »K*1 , 12) 

80   PRINT  9005, (Z  NU  (K)  f  Kal ,  16)  $  PRINT  9006»  ((3R  pR   <K)«K»ltl6) 

85   PRINT  9009,  (TMP(K) ,K=1, 32)  SPRINT  9010,  (S (K) ,K»1 ,32) 

9000  FORMAT (5X6HPSTART,6X4HDIAM,4X6HFACT0R»3X7HAVGFREQ, 
1  8X2HPM,9X1HN,7X3HNJJ.7X3HNFR»8X2HJD/ 

1  1XF10.2«F10.4,F10.6,F10.4,F10.0,2I10,F10*2,U0) 

9001  FORMAT (//1X7HTMR  VIS,2F7.2, 14F8.2) 

9002  F0RMATQX9HT  CON  GAS»F5.2»F7.2»14F8,2> 

9003  FORMAT (1X9HVIS  (E06) . 1 (6PF5.2) , 1 (6pF7.2) » 14 (6pF8.2) ) 

9004  FORMAT  (1X6HCONGAS,F8.0»F7.0»14F8.0) 

9005  F0RMAT(1X6HZ  NU   ,F7.3» 15F8.3) 
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9006  F0RMATUX6HGR  PR  ,F7.3» 15F8.3) 

9009  FORMAT (1X6HTMP    ,F7.3» 15F8.3) 

9010  F0RMATUX6HS      ,F7.3»  15F8.3) 
99999  END  READ  DATA 

C 
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19.     Appendix  F.     Subroutine  ELEMENT  for  Amplitude 


SUBROUTINE  ELEMENT 
CONSTANT  LENGTH  OF  EACH  ELEMENT  TYPE 
CALCULATION  OF  ELEMENT  LENGTHS  AND  TEMPERATURES 

C0MM0N/l/PSTART»0IAM»FRFtFACTOR»AVGFREQ,PM,N.NJJ,PI,AF,XKfCPGtM, 

1  TMP(32)»S(32)»TMRVIS(16) . VIS <16) tTCONGAS < 12) tC0NGAS(12) ,ZNU(16) . 

1  GRPR<16>iCVG»RH0DEWtTDEW»R,jD 

C0MM0N/2/TEMP(97) ,XX(97> »RH0(96) »XMID(96> »GM(96> tTM6(96> 
TYPE  REAL  NFR 

5  TEMP(1)=TMP(1)  S  XX(1)«S<1)  SNP1«N+1  $MP1=M*1 

10  0XXe(S(32)-S(l) )/N 

15  DO  40000  1*1, N   $XX(I*1)«XX(I)*D  XX 

20  DO  40001  K*2,32  $IF <S (K) .GE.XX (1*1 ) ) Go  TO  33 

40001  CONTINUE  $K*32 

33  TEMP(I*l>sTMP<K)*<TMP<K)-TMP<K-l)>/<S<K)-S<K-l))»<XX<I*l)-s<K>> 

35  TMG(I)«(TEMP(I)*TEMP(I*1))»0»5  $RHO ( I ) aPSTART/R/TMG ( I ) 

40  X  MID(I)«(XX(I)*XX(I*D)*0.5  $GM < I ) *RHo ( I ) #AF»DXX 

40000  CONTINUE 

45  DO  40002  I*NP1 »M$TEMP ( I*J ) «TMG ( I ) *T    DEW 

50  XX(X*1)"XX(I)*D  XX  $XMID{I)«(XX(I)*XX(I*1))»0.5  $RHO(I)=RHO  DEW 

55  GM(I)sRHO(I)*AF«D  XX 

40002  CONTINUE 
60  PRINT  9310 

9310  FORMAT (/) 

65  PRINT  8038* (TEMP(K),Kel.NPl) 

8038  FORMATOX  4HTEMR,  10F8.2»7F7.2) 
70  PRINT  8039, (XX <K) »K=1 »NP1 ) 

8039  FORMATOX  2HXX,  10F8.2»7F7.2) 
75  PRINT  9206,(X  MID (K) »K*1 ,N) 

9206  F0RMAT(4X  5HX  MID, 10F8.2»6F7.2) 
80  PRINT  9205* (TMG(K),K»ltN) 

9205  FORMAT (6X  3HTMG, 10F8.2»6F7,2) 

85  PRINT  9000»(GM(K) ,K«1»N) 

9000  FORMAT (7X2HGM,8E12. 4) 

90  PRINT  9207i(RHO(K),K=l»N) 

9207  FORMAT (6X3HRH0,8E12. 4) 

9208  F0rMAt(8X1HC,8E12.4) 

95  PRINT  9310  $PRINT8038» (TEMP (K) , KrNPl , MP1 ) 

100  PRINT  8039»(XX(K) ,K«NPl,MPl)  $  PRINT  9206» <X  MID (K) ,K=NP1 ,M) 

105  PRINT  9205* (TMG(K) ,K«NP1»M)  $  PRINt9000» <GM (K) »K=NP1 ,M) 

110  PRINT  9207t(RHO(K)fK"NPl.M) 

99999  END  ELEMENT 

C 
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20.  Appendix  G.  Subroutine  HEAT  TR 


SUBROUTINE  HEAT  TRUTESTFRF) 

COMMON/ 1/PSTART»DI AM »FKF, FACTOR, AVGFREQ,PM,N,NJJ, PI ,  AF,XK,CPG»M, 
1  TMp(32)fS(32)»TMRVlS(16) t VIS (16) »TC0NGAS(12) » CONGAS < 12) ,ZNU < 16) , 
1  GRPR<16) »CVG,RH00EW»TDEW,R,JD 

C0MM0N/2/TEMP(97) ,XX(97) »RH0(96) »XMID<96)  *GM(96) ,TMG(96) 

C0MM0N/3/G(2) ,P<2) ,T(   360) ,PRES (2) ,X (2,360) ,XAVG(2) , EKE (360). 
1  VOL(    2),R0(    2) ,FR(96) ,ZKE(96) tQ(96) ,W0RK(96) ,IFRF, 
1  VEL,R0E,THETA,AREA,DQ,DT,TWI,TAVG,REYN0,HHT1,HHT2,HHT 

C0MM0N/4/NFR 

TYPE  REAL  NFR 

1  IFdTESTFRF.EQ.DGO  TO  386 

2  DO  40014  K=2,12  $IF(T  CON  GAS(K).GE.T  AVG  ) GO  TO  3 
40014   CONTINUE  $K=12 

3  CON  G=C0N  GAS(K)*(CON  GAS(K)-CON  GAS(K-1  ))/(T  CON  GAS(K>- 
1  T  CON  GAS(K-1  ))<MT  AVG  -J  CON  GaS(K)) 

4  DO  40013  K=2,16  $IF<TMR  VlS(K).GE.T  *VG  ) GO  TO  5 
40013   CONTINUE  $K=16 

5  VISC=VIS(K)+(VIS(K)-VIS(K-1  ))/(TMR  VlS(K)-TMR  VIS(K-1  ))* 
1  (T  AVG  -TMR  VIS(K) ) 

6  IF(VISC.EQ.0.)GO  TO  400 
8   REYNO=ABSF(DIAM»VEL*ROE/VISC) 

10  IF(REYNO.LT.0.00l)REYNO=0.001 

20   CONTINUE 

22  IF(CONG)4Q0,400»25 

25  PRANDTL=CPG     *VISC/C0N  G 

26  IF(PRANDTL)  400*400*30 

30  IF(TAVG)400, 400,31 

31  ZNUSSELTso,023*REYNO#*0.8»PRANDTL**0.4 

32  HHyl=      ZNussELj*CON   G/DIAM#FACtOr 

33  THETA=TWI-T    AVG      SIFUHETA)    35,385,35 

35   GP=CP6*ROE**2*980.66/T  AVG*DIAM**3/VISC/C0NG*ABSF (THETA) 
37  IF(GP)  410,40,40 

40  IF(GP.EQ.0.)ZNUS=0. 

41  IF(GP.EQ.O.)GO  TO  380 

42  GP  L=LoGF(Gp)»0. 43429 

43  DO  40016  K=2,16  $IF(GR  PR(K).GE.GP  L)GO  TO  370 
4OOI6   CONTINUE  $K=16 

370  ZNUL=ZNU(K)*(ZNU(K)-ZNU(K-1  ))/ 

1  <GR  PR(K)-GR  PR(K-1  ))*(GP  L-GR  PR(K)) 

375  ZNUS=10.«»ZNUL 

380  HHT2=ZNUS«C0N  G/DI AM*FACTOR  $HHT    =HHT1    +  HHT2 

385  DQ=HHT»AREA*THETA*DT 

386  CONTINUE  $IF(      REYNO  .LT. 0.001 ) REYNO»0.001 
2386  IFdFRF.EQ.DGO  TO  99999 

387  IF(REYNO.LT.1000)FRF=16. /REYNO  »NFR 

388  IF(REYNO.LT.1000)GO  TO  99999 

389  FRF=0.01«NFR 

390  FRF1=FRF 

392  FRF=(l./(3.2*0.43429*LOGF(REYNO*SQRTF(FRFD >*1.2> ) **2*NFR 

393  IF(ABSF(  (FRF1-FRF)/FRF)  .GT.0.000DGO  To  390 
395  GO  TO  99999 
400  PRINT  9000.CON  G, VISC,PRANDTL,TWI ,T  AVq 

9000  FORMAT (//I6X5HCON  G» I6X4HVISC, 13X7HPRANDTL, 17X3HTWI , I5X5HT  AVG/ 
1  1X5E20.5) 

410   PRINT  9001»CON  G,  VISCGP,  TWI  ,T  AVG 

9001  FORMAT (//16X5HC0N  G, 16X4HVISC, 18X2HGP      , 17X3HTWI , 15X5HT  AVG/ 
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1  1X5E20.5) 

420  DT=0. 

99999  END  HEAT  TR 
C 
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21.  Appendix  H.  Subroutine  AMPLITUD 


COMPUT 


1 

2 

3 

40009 

9004 

9005 

4 

5 

40004 

6 

10 

15 

25 

30 

31 

32 

35 

40008 

40 

45 

50 

55 

60 

61 

65 

70 

80 

85 

90 

95 

100 

105 

110 

120 

240 

241 

245 

250 

255 

260 

265 

270 

275 


SUBROUTINE  AMPLITUD 
ATION  ASSUMES  PRESSURE  VARIATION  IS  A  COSINE  FUNCTION 

DIMENSION  PX(96) ♦KOUNTMAX (96) »JJMAX(96) 

COMMON/1 /PSTART.DI AM* FRF ♦ FACTOR ♦AVGFREQfPM,N,Nj J, PI »AF«XK.CPG,M, 

TMp(32)»S(32)  1TMRVISU6)  .VIS(16)  tTC0NSAS(12)  »CONGAS(12)  ,ZNU(16)  , 

GRPRU6)  »CVG»RHODEW»TDEW»R,jD 

C0MM0N/2/TEMP(97) ,XX(97) »RH0<96) »XMI0(96>  »GM<96) »TMG(96) 

C0MM0N/3/G(2) »P(2)*T(   360) ,PRES(2> »X (2,360) t XAVG (2) *EKE  (360) . 

VOL(    2),R0(    2)»FR(96) ,ZKE(96)»Q(96) ,W0RK(96)»IFRFf 

VEL, ROE* THETA, AREA »DQ,UT,TWI»TAVG,REYNO,HHTl»HHT2»HHT 

TYPE  REAL  NFR 
IF(FRF.NE.0.)IFRF=1  $IF (FRF.EQ.0. ) IFRF=0 

PRINT  9004 

DO  40009  K*l»7  $REYNO=10.*»K  $ITESTFRF*1  SCALL  HEAT  TR(ITESTFRF) 

PRINT  9005»REYNO,FRF 

CONTINUE 

FORMAT (//10X,5X5HREYNo»7X3HFRF) 

FORMAT(lOXf2El0.3) 

PRINT  9006 

PRINT  9006 
DO  40004  J*l,360  $X (2» J) «XX d )  $EKE(J)aO.O 

CONTINUE 

NJJL2*NJJ-2 

DT»JD/360./AVGFREQ   $NL1«N-1  $  NP1=N*1     $  SUM  W0RK=SUM  Q=0. 

SUMFRsGlD2aSUMEKEaO.  $ROlD2eRHODEW 

MP1*M*1  SXFaXX(N*f)  $C8sPl«DlAM  $C9=»  (XK-1 .  ) /XK   $XL=XF-XX(1> 
DO  40000  Ial.M 

JJMAX(I)»0 

G(2)*GM(I) 
DO   40008   K-1O60    $ X (1 »K) -X (2»K)       $T(K)*TMG(I) 

CONTINUE 

IQ=0 

V0L(2)»GM(I)/RH0(I) 

X(2*360)»X(1*360)*VOL(2)/AF 

R0(2)=RHo(I)  $XAVG(2)=(X(l,360)+X(2,360))*0.5$Q(I)=0. 
DO  40002  JJsl,NJJ  $IF(IQ.GT.2)G0  TO  40002 
IF(JJ.GT.JJMAXd)  )  JJMAX(I)«JJ 

QI»Q(I)  $KOUNTMAX(I)=0  JWORKl«WORK ( I ) 

FR(I)=    ZKE(I)=Q(I)=WORK(I)«0, 

HO  40003  J=JD,360,JD  $R2«J»Pl/180«$SlN2xSINF <R2) 

IF(J.EQ.JD) JLJD=360  $IF ( J.GT. JD) JLJD=J-JD 

RO(   1>*R0(   2)  *VOL(   l)*VoL(   2) 

XAVG(   l)aXAVG(   2)   $R0lDl»R0ID2 

PRES(   1)=PRES(   2)  $G(   1 )  sG  (   2)  $  P   <   1 > "P   <   2) 

KOUNTsO      $GID1=GID2 
IF(XUtJ)  ,GE.XF.AND.X(1»JLJD)  .GE.XF)GO  TO  690 

DQ=AQ       so. 

CONTINUE  $KOUNT=KOUNT*l  $Dq1»Aq 
IF (KOUNT.GT.KOUNTMAX ( I ) ) KOUNTMAX ( I ) aKOUNT 

X(   2, J  )sX(ifJ  )*VOL<   2J/AF 

XAVG(   2)«(X(   2,J  )*X(1,J  ))*0.5 

vEL=(XAVG(   2)-XAvG<   1 ) ) /Dj 

TAVG=(T(JLJD)*T(J) )*0.5 
IF (XAVG (2) ,LT.XF)PX(I)=PM*C0SF(PI/2.»(XAVG(2)-XX(1))/XL) 
IF(XAVG(   2).GE.XF)PX(I)«0. 

P(   2)sPX(I)#SIN2  $PRES(   2)-P<   2)*PSTART 
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280   PRESS=(PRES(   1>*PRES(   2) ) »0.5  $ROE= (RO (   2)*RO(   1>)»0.5 

300   IF<x(2»JLJD).GT.xF.0R.X(2»J).6T.xF)G0  TO  500 

305   V0LM=6M(I)/R0E 

400   AREAsC8*(X(2»J)-X(l»J)*X<2,jLJD)-X(l,JLJD) )«0.5 

404  G(2)aGM(I) 

405  AVGXs(XAVG(   1)*  XAVG(   2))  *0.5 

410    DO  40005   K«»2»NP1  $IF (XX (K) .GE.AVGX) Go  TO  420 

40005  CONTINUE  $K«NP1 

420   TWI=TEMP(K)* (TEMP (K) -TEMP (K-l) ) / (XX <K) -XX (K-l ) ) » (AVGX-XX <K) ) 

425       ITESTFRFbO  SCALL  HEAT  TRdTEST  FRF) 

430  IF(DT.EQ.0.)GO  TO  99999 

435   0  FR=ABSF(FR  F/2.  »AREA*R0E*VEL*»3»DT) 

440   DWK     «C9»DQ-V0LM/XK*(P(   2)-P(   1 ) ) -OFR/XK 

455   QWaDQ-DWK        $AQSDQ 

460   T(J)=T(JLJD)*QW/GM(I)/CVG 

465   RO(   2)"PRES(   2)/R/T<   J) 

470   VOL(   2)«GM(I)/RO(   2) 

475  IF(K0UNT.LT.2)G0  To  240 

480  IF(KOUNT.GT.10.)GO  TO  700 

485  IF(AQ.EQ.0.)AQ*0.00000l  $  IF  (ABSF ( (AQ-DQ! ) /AQ) .GT.O.OOOl  ) GO  T0240 

490   GO  TO  700 

500  IF(X(   2»J  ) .GE.XF)C2=XF  $IF(X(   2»J  )  .LT.XF) C2«X <   2. J  ) 

510  IF(X(1   tJ  ).GE.XF)W2*XF  $IF(X(1   ij  ) .LT.XF) W2=X ( 1   .J  ) 

520   lF(X(2tJLJ0).GE.XF)Cl=XF  $IF (X (2f JLJD) .LT.XF) Cl=X <2» JLJD) 

530  IF(X(1,JLJD).GE.XF)W1=XF  $IF (X (1 , JLJD) •LT»XF) W1=X ( 1 , JLJD) 

550  IF(VEL.GT.0.0)EKE(J)»0.5*(VEL*AF»ROE»DT)*VEL*VEL 

551  IF(VEL.LE.0.0)EKE(J>*0.0 

600   AREAls(Cl-Wl)#C8  $  AREA2« (C2-W2) #C8  $AREA= (AREAl*AREA2) #0.5 

605   AVGX=0.25»(W2*W1+C2*C1) 

610    DO  40006  KR-1.NL1  $K=N  -KR  $IF  (XX (K) .LE. AVGX) GO  TO  620 

40006  CONTINUE 

620   TWI=TEMP(K)*(TEMP(K)-TEMP(K*1))/(XX(K)-XX(K*1))*(AVGX-XX(K)) 
625       ITEST  FRF*0  SCALL  HEAT  TR(ITESTFRF) 
630  IF(DT.EQ.0.0)GO  TO  99999 

635  0  FR=ABSF(FR  F/2.  *AREA*R0E*VEL*«3*DT> 

639  V0LM=(G(1)*G(2) )»0.5/ROE 

640  DWK     aC9*DQ-V0LM/XK«(P(   2)-P(   l))-OFR/XK 
655  QW=DQ-DWK        $AQ*Dq 

670   V0LIP2=AREA2»DIAM/4.  $G (2) «RO (2) *V0LIP2  $GID2»GM ( I ) -G (2) 

674  IF(GID2.GT.GIDl.AND.GID2.GT.0t) 

1  R0ID2=(R0ID1«GID1*R0(1)*(GID2-GID1))/GID2 

675  IF(GID2.LE.GIDl)R0ID2=R0IDl 

676  V0LID2SGID2/R0ID2    SVOL <2) aVOLIP2*VOLlD2 

680  IF(G(   2).LE.0.)T(   J)»T  DEW  $  IF(G(   2).LE.0.)GO  TO  687 

681  IF(G(2)-G(1>)685.685»683 

683  T<J)=(T(JLJD)«G(1)*TDEW»(G(2)-G(1) ) ) /G (2) *QW/CVG/G (2) 

684  GO  TO  687 

685  T(J)=T(JLJD)*QW/CVG/G(2) 

687  CONTINUE 

688  RO(   2)*PRES(   2)/R/T(   J)  $IF (K0UNT.LT.2  ) GO  TO  240 

2688  IF(KOUNT.GT. 10)60  TO  700 

689  IF(Aq.Eq.O.)Aq«0. 000001  $  IF (ABsF ( (Aq-Dq1 ) /AQ) .GT. 0.0001  > GO  T0240 

2689  GO  TO  700 

690  DQ=DQl=AQsDWK=QW=D  FRaQ,  $REYNO«.001  $HHTlaHHT2aHHTaO. 

691  ITEST  FRFal  SCALLHEAT  TRdTEST  FRF) 

692  G(   2)a0.  $T(   J)aT  DEW  $r0(   2)«R0ID2    $VOL   (2) aQM ( I) /RO (   2) 
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693  X<   2»J  )=X(1,J  )*VOL<   2)/AF  $XAV6<   2)=<X<2   ,J  )*X<1,J  ))*0.5 

694  VEL=(XAV6(   2)-XAVG(   1))/DT 

695  PRES(   2)=PSTART 

696  TWI=0.  $K0UNT=K0UNT*1 

697  IF(KOUNT.GT.KOUNTMAXd)  )  KOUNTMAX  ( I )  =KOUNT 

700  CONTINUE 

701  WORK(I)=WoRK<I)*DWK  $  U(I)bQ(I)*DQ  $FR < I ) =FR (I ) *D  FR 

703  IF(VEL.GT.o.)DKE=0.5*(VEL»AF«ROE*DT)*VEL«VEL 

704  IF(VEL.LE.0.)OKE=0.  $ZKt < I ) sZKE  (I ) *DKE 

705  IF(IQ.LT.1)G0  To  40003 

711  IF(J.EQ.90)PRINT  9001 
9001   F0RMAT(///3XlHI,3XlHJ,iX2HJJ»lX5HK0UNT,   5X6H   P<2)»4X7H   R0(2)» 

1  5X6H   T(J) »5X3HTWI,8X3HVEL.9X2H0Q»8X3HDWK,8X3HDFR»3X8H   X(2,J), 
1  1X8H   V0L(2)»5X6H   G (2> /63X ,6X5HREYNO» 7X4HHHT1 . 
1  7X4HHHT2.8X3HHHT.8X3HFRF) 

712  IF<J.EQ.90.Or.J.EQ.180.Or.J,EQ.270.Or.J.EQ.360> 
1  PRINT  9000tI.J»JJfKOUNT»P   (   2)»RO(   2)»T(   J) ,TWI ♦ VEL»DQ,DWK. 
1  DFR,X(2»J)»V0L(2) ,G(2) .  REYNO.HHT1. 
1  HHT2tHHT,FRF 

9000   F0RMAT(lXI3»I4.l3,l6,Fll.l,E11.4,Fii,5,F8.3,4EH.3iFll.5»F9.4, 
1  E11.4/63Xt6Ell.4) 
720  IF(J.EQ.360)PRINT  9003.U < I ) , WORK  ( I ) » ZKE  ( I )  »FR(I) 
1  ,QLWOQ»DELQIOQI»OELWKOWK 
9003   F0RMAT(/12X4HQ(I) ♦  8X7HWORKU).   9X6HZKE < I ) » 10X5HFR < I ) 
1  ,10X5HQLWOQ»7X8HOELQIOQI»7X8HDELWKOWK 
1  /1X4F15.3.3F15.5/) 
40003   CONTINUE 

723  IF(JJ.LT.2)G0  TO  40002 

724  QI_WOQ=DELQIOQI=DELWKOWK  =  0.0 

725  IF(Q(I) .Eq.O.O.AND.WoHK(I) ,Eq.0.0)Go    To    750 

726  IF(Q(I) .EQ.0.G)Q< I) =0.000001  S IF  (WORK (I) .EQ. 0.0) WORK < I ) =0.000001 
730   QLWOQ=ABSF( (Q ( I > -WORK ( I ) )/Q(I) )  $DELQlOQl=ABSF ( (Ql-Q ( I ) ) /Q ( I ) ) 
735   DELWKOWKsABSF< (WOrKI-WOrK ( I ) )/WOrK(I) ) 

738  IF(QLWOQ.LT. 0.0001. AND. DELQIOQI.LT. 0.0001)60  TO  750 

739  IF(QLWOQ.LT.0.000l.AND.DELWKOWK,LT.0,000l)G0  TO  750 

740  IF<DELQlOQI.LT. 0.0001. AND. OELWK0WK.LT.O.OOODG0    To    750 
742    IF(JJ.GE.MJJL2)G0    TO    750 
745      GO    TO    40002 
750       IQ=IQ*1 

40002   CONTINUE 

40000  CONTINUE 

756  PRINT  9006 
9006   FORMaT(IHI) 

757  PRINT  9501 
9501   F0RMAT(1X,4X1HI.1X8HK0UNTMAX»1X5HJJMAX, 

1  5X5HPX(I) ,3X7HWORK(I) » 3X7HSUMW0RK.6X4H0 ( I ) . 

1  6X4HSUMQ.5X5HFR(I) ♦5X5HSUMFR*4X6HZKE < I ) ) 

760  DO  40001  1=1, M  $SUMq=SUMq*q < I )  $SUM  WORK=SUMWORK*WORK ( I ) 

761  IF  (SUMWORK.EQ.0.)SUMWORK  =  0. 000001 
765   SUMFR=SUMFR+FR(I) 
775   PRINT  9500»I»KOUNTMAX(I) ♦JJMAX(I) »PX(I) 

1  ♦W0RK(I),SUMW0RK,Q(I) »SUMq,FR(I) »SUMFR, ZKE ( I ) 
9500   FORMAT(1X,I5,I9.I6,8F10.0) 

40001  CONTINUE 

780    DO  40007  J=l,360  $SUMEKE=SUMEKE*EKE < J) 
9310   FORMAT(//) 
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40007   CONTINUE 

785   EKEOSW=SUMEKF7SUMWORK  $ZKENOSw=ZKE (N) /SUMwORK 
787   ZKEMOSW=ZKE(M) /SUMWORK 
PRINT  9008 

9008  FORMAT (///1X»9X6HSUMEKE»9X6HEKE0SW,8X7HZKEN0SW,8X7HZKEM0SW) 
PRINT  9007»SUMEKF.»EKEOSW»ZKENOSW»ZKEMOSW 

9007   fORMAT(1X,4F15.4) 
PRINT  9006 

Do40010  J=JD«360»JD 
PRINT  9009fJ,EKE(J) 

9009  FORMAT (1X,6HEKE (J) ,  I3»EU,4) 
40010   CONTINUE 

99999   END  AMPLITUD 
C 
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22.  Appendix  I.  Subroutine  READ  DATA  for  Frequency 


SUBROUTINE  READ  DATA 
COMBINATION  OF  STATEMENTS  FOR  FREQUENCY  TYPE  READ  DATA 

COMMON/l/N»TMP(32) ,S(3?) ,RH0DEW» TDEW.M.XK.R. TEMP (97) .XX(97) *PIi 
1  C(96) «TMG(96) »XMI0(9fc) .RH0<96) .GM (9ft) . A (9ft) t Z (96) .AFtPSTARTt 
1  DXXOVC 

READ  8001tPSTART»DlAM,M 
8001   FORMAT(2F10.C»1 110) 

RFAD  8000*  (TMP(I) tl  =  i»32)  $READ  8000» <S  ( I ) » 1  =  1 .32) 
8000   FORMAT(16F5.0) 
CONVERSION  OF  INPUT  DATA  TO  UNITS  USED  IN  CALCUl  ATIONS 
COMPUTATION  USES  CENTIMETERS  RUT  INPUT  TS  IN  INCHES 

DO  40012  K=1.32  $S(K)=S(K)*2.54 
40012   CONTTNUE  $DI AMbDI AM*2.54 

WT=4.003SRU=8.317E07$R=RU/WT$PI=3.1415926S36$AF=PI/4,«DIAm**2 
XK=5./3.  $TDEW=4.3  $m=n*N/2  $RHODEWsPSTART/R/TDEW 
COMPUTATIONAL  VALUES  ARE  PRINTED 
PRINT  9000»PSTART.DIAM,N 
9000   FORMAT (1X13HPSTART»DI AM. N»F8.0»F8.4» 18) 

PRINT  QQo9« <TMP(K) »K=i ,3?)  SPRINT  90l0 * <S (K) .K  =  l ♦ 32) 

9009  F0RMAT(1X5HTMP   16F8.2) 

9010  F0RMAT(1X5H   S   16F8.2) 
99999   END  REAO  DATA 

C 
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23.     Appendix  J.     Subroutine  ELEMENT  for  Frequency 

subroutine:  element 

common/1 /nttmp (32) »s(32) .rhodew* tdew»m, xk »r» temp (97) .xx (97) » pi ♦ 

1  C(96)  »TMG(96)  »XMID(9ft)  »RH0(96)  »GM <9ft ) » A (9ft ) » Z (96)  .AF.PSTART. 

1  DXXOvC     

CALCULATION  OF  ELEMENT  LENGTHS  AND  TEMPERATURES 
10   TEMP(1)=TMP(1)  $  XX(1)=S(1)  $NP1=N*1 
?0      0  XX  OV  C»(S<ip)-S(l))/N/SQRTF(XK*R*TMp(16)) 

24  00    40002    I=1»N 

25  C(I)»SORTF(XK*p*TFMP(I)) 
?7      CONTINUE 

28  XX(I*1)=XX(I)*D    XX   OV   c*c<n 

29  00  40000  K=2.32  5  IF (S <K) .GE.XX ( I *1 ) > GO  TO  33 
40000   CONTINUE  *  K=3? 

33   TFMP(I*1)=TMP(k)*<TMP<k>-TMP(k-  1  > > / (S (K) -S (k-  1 > > • <XX <I*1 ) -S («>> 
17   TMG<T)=(TEMP<I)*TEMP(I+1) )/2.  $  CCCsSQRTF  (XK*R*TmG  ( I )  ) 

38  IF(ABqF(CCC-C<I) )/C(I)-.  001)43*43*39 

39  C(I)aCCC  $  GO  TO  27 
43   C(I)=CCC 

4000?   CONTINUE 

47  IF(AB«!F(XX(N*1)-S(32)  )  /S  <32)  -.0001 )  50  »5o  .49 

49  O  XX  OV  C=D  XX  OV  C*(1.*\.*<(S(3?)-XX(N*1))/XX(N+1) ))  SGOTO  24 

50  DO  40017  I  =  1«N  $X  MIO(I)  =  (XX(I)*XX(I  +  l)  )/2, 

SS   RHO(T)=P  START/R/TMG(I)  SZ < I ) =C < I ) *RH0 ( I )  $A ( I ) =1 ,/Z <  I ) 
ftO   GM(I)=RHO(I)*AF    * (XX ( I +  1 ) -XX  ( I )  ) 
40017   CONTINUE 

ss  print  9310 

9310   FORMAT*/) 

<?0      PRINT    9038.  <TEMP(K)  iK=liNPD 

8038  format (ix  4HTEmp«ioP8.2»7F7,2) 

9*  PRINT  8039«  (XX(K)  »K*] »NPi)   ! 

8019  FORMATOX  2HXX.  10F8.2.7F7.2) 

lnft  PRINT  R2o6*<x  MID(K) »K=1.N) 

92n6  FORMAK4X  5HX  MID*  10F8.  2,  6F7.2) 

loS  PRINT  9205» (TMG(K) iK=i»N) 

92oS  FORMAT (6X  3HTMO. i0F8.2,6F7.2) 

110  PRINT  R000» <GM(K) .K=] *N) 

9000  F"ORMAT(7X2HGM»8El2.4) 

120  PRINT  9207» (RHO(K) iK=l»N)  $  PRINT  92o8» <C (K) • Kel .N) 

9207  FnRMAT(6X3HRH0»8El2.4) 

9208  F0PMAT(8X1HC»8F12.4) 
EVD    FLFMENT 
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24.  Appendix  K.   Program  FREQENCY 


PPOGpAM  FREQENCY 
C  PROGRAM  OF  JUI  Y  1967   WITH  OMISSIONS  AND  PRINTOUT  DELFTIOnS  OF  OCT  6  1967 

COMMON/ 1/N.TMP (32) tS(32) ,RHODEW»TDEW»M, XK.R. TFMP (97) »XX(97) »PI» 
1  C(96)  «TMG(96) »XMI0(9ft) ,RH0(96) .Gm (96) ♦ A (96) .Z (96) .AF.PSTART. 
1  DXXOvC 

COMMON/2/P(96,]02) »U<96»102) .PMAX(96*12) .UMAX (96. 1 2) .UMlN (96. 12) . 
1  TMIN(96»12)  .TMAX(96H2)  ♦dmIN(96«1?) 
l.ZlNTFRVLd?)  »PFRI0D(12)  »FREQ(12) 

5  call  read  data 

10  call  element 

15  Print  9310 

?0  kounfs-100  $  sambda^ooq.  $gambda=0 . 1*sambda 

?5  n.j=1^2  snj  less  1=nj-1  $  nj  less  2=nj-2$nl.1=n-1  $n  plus  1=n*1 

30   DO  4o000J=l*2  $DO  40000  I=1*N  $  P d » j) =-SAMBDA  $  U(I.J)=0. 

40000  CONTINUE  5  DO  40001  JKsl.10  $  DO  40001  I=1»N 

^«;     p  max<i»jk)=u  max(i»jk)=ii  min  ( i  » jk  )  =0  .$t  min(i.jk)*t  max(i.jk)=1. 
40     ifuk.fq.1jp  min(i.jk)=-samboa  $if  ( jk  .gt,  1  )  p  min  ( i .  jk)  =sambda 

40001  continue  $  jk=i   $  max  coiinT=o 
calculation  starts 

50  u(n.?)=(-1)*sambda*a(n)  %   p<n.2)=0. 
60  kount=koumt*l00  s  sum  frfq=0. 

65  do  141  j=3»nj  *j2=j+k0unt  $if  (  < j/2) »2-j) 101 » 121 » 101 
computation  of  odd  time  pressure  and  velocity 

101  DO  115  »I=2»N,2 

105  P(L  D=P(I«J-2)+2.<1'(U(I-i.J-l)-U(I»J-l5  )/(A(I-i)*A(IJ) 

107  P(T-i *J)=P(I»J) 

111  U<i»J)=-U(I»J-2)*2.<MlJ<I-T  •  J-l)  *A  (I)  ♦U(I»  J-1  )*A  ( i-l)  ) 

1  /(A(T-1  )*A(I) ) 

113  Ud-i .J)=U(I.J) 

21J4  IF(P(I.J) .GT.P  MAX(I»JK))T  MAXdtJKjsJ? 

2115  iF(Pd.J)  .GT.P    MAX(I,JK))P    MAX  ( I ,  JK)  =P  ( I .  J) 

2116  lF(Ud,J)  .GT.U    MAXd.JKMU    MAX  ( I , JK) =U ( I , j) 

2117  iF(Pd.J)  .I.T.P    MIM(I,JK))T    MIN(I.JK)=J2 
2iTfl  iF(Pd.J)  .LT.P    MIN(I»JK))P    MIN(I.JK)=P(I.J) 

2119  IF(Ud.J)  .LT.U  MlMd.JKJlU  MTN  d  .  JK)  sU  ( I .  J) 

2120  iF(Pd-l.J)  .GT.P  MAX(  I-]tJK))T  MAX  ( I~l  .  JK  )  =J2 

2121  iF(Pd-l.J)  .GT.P  MAXd   -1»JK))P  MAX  (  1-1  ♦  JK)  =P  ( 1-1 .  J) 

2122  lF(UfI-l  .J). GT.U  MAX  ( 1-1  » JK ) ) U  MAX ( I"l . JK ) =U ( 1-1 . J) 

2123  iF(Pd-l.J)  .LT.P  MINd-l.JK))T  MIN  ( 1-1 .  JK)  =J2 
21?4  IF(P(I-1.J)  .LT.P  MlNd-l.JK)  )P  MIN  <I-1  .  JK)  =P  (1-1 ,  J) 
2125  IF(Ud-l.J)  .LT.U  MIN  H-l  •  JK)  JU  MIN  ( t-1 .  JK)  =U  ( 1-1 ,  J) 
21P6  IF(I.GT.?)G0  TO  115 
2127  IF(P(I.J) .GT.GAmBOA)mAX  COUNT=MAx  COUNT+1 

212B  IF(MAX  COUNT. Gf.N. AND. P(1»J) .LT.O.     )JK=JK*1 

2129  IF(M*X  COUNT. GT.N. AND. P(l »J) .LT.O.     )  MAX  CDUNT=0 

115  CONTINUE 

117  GO  TO  141 
COMPUTATION  OF  EVEN  TIME  PRESSURE  AND  VELOCITY 

121  Ud.|)=0.  SP<ltJ)=Pd»J-2)-2.*U(2»J-i)/A(l) 

3122  IF(P(1,J) .GT.P  MAX(1,JK))T  MAX(1,JK)=J2 

31P3  IF(Pd.J)  .GT.P  MAX(1,JK))P  MAX ( 1 , JK) =P d » J) 

3124  IF(U(1  »J) .GT.U  MAX(1,JK))U  MAX ( 1 , JK) =U (1 > J) 

31?5  iF(PlliJ) .LT.P  MIN(ltJK))T  MIN(1,JK)=J2 

31?6  IF(P(1.J) .LT.P  MIM(1,JK))P  MTN(1.JK)=P(1.J) 

3l?7  TF(Ud,J)  .LT.U  MIN«1,JK))U  MTN  ( 1  ♦  JK  )  =U  (  j  ,  j) 

31?B  IF (P(l . J) .GT.GAMBOA)MAX  COUNT=MAX  COUNT*) 
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3129  IF(MaX  COUNT. GT.N, AND. P<1 »J) .LT.O.     )JK=JK+1 

3130  IF(MAX  COUNT. GT. M. AND. P(nJ)  .LT.O.     )  MAX  COUNT=0 
1?5  00  134  1=3.  N|_]  »2 

1?7  P(T»J)=P(I»J-2)*2.»(  U(I-ltJ-l)-U(I,J-l)  )/(A(I-l)*A(I)  ) 

129  P{I-l.J>=pn»J) 

131  U(T.J)=-U(I»J-2)*2.*(U(I-i  tJ-i)*A{I)*U(I,j-i)*A<I-i>) 
!  /(A(T-i  )*A(D  ) 

132  U(I-HJ)=Ud'J) 

4133  IF(P(I.J)  .GT.P  MAX<I,JK))T  MAX(I.JK)=J2 

4P4  IF(P(T«J)  .GT.P  MAX(I.JK))P  MAX ( I , JK) =P ( I , J) 

4135  IF(UU.J)  .GT.U  MAXlltJKMU  MAX  ( I ,  JK  )  oU  ( I  ,  j) 

4136  IF(P(I»J).LT.P  MZN(X«JK))T  MIN(I,JK)=J2 

4137  IF(P«T*J) .LT-P  MlN(ItJK))P  MTN < I . JK ) =P ( I , J) 
41T8  IF<U(I,J) .LT.U  MIN(I»JK))U  MIN ( I » JK) =U ( I . J) 

4139  IF(P(I-1.J) .GT.P  MAX(  I-],JK))T  MAX  ( 1-1 ♦ JK) =J2 

4140  IF(P(I-1.J) .GT.P  MAX(I   -1 . JK) ) P  MAX (  1-1 , JK) =p ( 1-1 , J) 

4141  lF(U(i-l  »J) .GT.U  MAX<I-1 tJK) )U  MAX ( 1-1 ♦ JK ) =U ( 1-1 »J) 

4142  IF(P(I-1»J) .LT.P  MIN(I-I.JK) )T  MIN < 1-1 . JK) =J2 

4143  IF(P(I-1'J).LT.P  MIN(I-1,JK)  )P  MIN  (  T-l  .  JK)  =P  ( 1-1 1  J) 

4144  lF(IJ(I-i»J>  .LT.U  MIN(I-I.JK)  )U  MTN ( T-l , JK ) =U < 1-1 » J) 

134  CONTINUE 

135  P(N.  |)rO.      $U(N»J)=U(N»J-1) *(P<N.J-1)-P(N.J) )*A(N) 

5136  IFfP(N.j)  .GT.P  MA)UN,JK>)T  MAX(N,JK)=J2 

5137  iF(PfN.J) .GT.P  MAX(N,JKMP  MAX <N, JK) =P (N, J) 
5l3fl  IF(U(N,J) .GT.U  MAX(N,JK))U  MAX (N, JK) =U (N, J) 
5n<)  lF(P(NtJ)  .LT.P  MIN(N,JK))T  MIN(N,JK)=J2 
5l4f)  IF(P<N,J)  .LT.P  MIN(N.JK))P  MIN (N, JK ) =P <N, J) 
5141  IF(U(N.J) .LT.U  MIN(N,JK))U  MIN <N, JK ) =U (N» J) 

141  CONTINUE 

150  IF(K0UNT.GT.  0  . AnO.KOUNT.LT. 1000 > GO  TO  226 

16n  00  2?5  J=1,NJ  LESS  2  $  J2=J+K0UNT 

220  IF(N.LF.16)PRIMT  9220 ♦ J2. <P ( I « J) . 1=1 »N) 

221  IF(n.lE.16)PRInT  9221 » (U(ItJ) .1=1. N) 

222  IF(N.GT.16)PRINT  9220 » J2, (P <I ♦ J) , 1=1 » 16) 
2?4  TF(N.GT.t6)PRTNT  9222» (P ( 1 • J) » I=17»N) 
2?3  lF(N.GT.l6)PFriNT  9221  »  <U  ( I  .  J)  .  \  =  \  »  16) 

5223  T.F(N.GT.l6)PRlMT  9221 ♦ <U < I • J) , I=\7 ,n) 

225  CONTINUE 

226  IFUK.LT.  11)60    TO    229 

6227  DO    4n003    J=l»9 

6228  ZINTFRVL(J)=T    MlN(l»J*l)-T    MIN(l.J) 
62?9  PFRIOO(J)=ZINTFR\/1  <J)*l)    XX    OV    C 

6230  FRFQ(j)=i./PERIOD(J) 

6231  SUM  fPFQ  =  SUM  FREO-t-FRFO(J) 
40003  CONTTNUE 

6232  AVG  FRFQ=SUM  FrEQ/9. 

6233  OMEGA=2.*Pr*AVG  EPEG 

6234  PRINT  93to 
P310  FORMAT (//) 

6235  DO  4«i005  J=l»l0 

6236  PRINT  9300»J  »(T  MIN ( I , J) ♦ 1=] ., 16   )  ,  <P  MIN ( I .J) ♦ 1=1 » 16   ), 
1  (IJ  MTN(I, J)  .1  =  1,16   )  ,  <T  MAXU,J)  »I  =  1»16   ). 

I  <P  MAX(I, J) ,1=1*16   ),<U  MAX(T,J) ,I=1»16   ) 

9300  FORMATM07H0MAXTMUM  AND  MTNIMUM  PRESSURES  AND  VELOCITIES  AND  CORRF 
1SPONDTNG  CALCULATION  TIMES  IN  ELEMENTS  1  THROUGH  16. 
11X12HTHTS  IS  THE  I2.9H  TH  CYCLE 
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11X5HT  MINIGFS.O./IXFHP  MIM16F8.0 » /1X5HU  MIM16F8.0./ 
1  1X5HT  MAX16F8.o«/lX5MP  MAX16F8  .0  •  /l  X.5HIJ  MAXlfFfi.O) 
400nS   CHNTTNUF 

6?39   PRINT  930?.  U»J  =  1»9)  ♦  <ZlNTEPVL<J)  »J=1 *9) ♦  (PFRIOD(J)  »J  =  1»9)  » 
1  tPRFn(.J)  »J=1»9)  ,AV<5  FRFQ.OMFGA    *D  XX  OV  C 

930?  FORmAT(131h6ZInTERVl  IS  TMF  nUmRER  OF  CA|CUiATIon  UnITS  In  OnF  CYC 
1LF.   PERIOD  IS  THE  TIME  LENGTH  OF  ONF  CYCLE  AMD  THE  RECIPROCAL  OF 
1THF  PFRTOn  /126H  I*;  THE  FREQUENCY.  THE  TOP  LINE  GIVES  THE  CYCLE  Nil 
1MRFR  FOLLOWING  THE  INITIAL  PULSE.  AVG  FQFO  IS  THE  AVERAGE  OF  ALL  F 
1RF0UEMCIES  / 

UXftHCYCl.E     9112,/lXRHZlwTEPVL  9Fl2. 0  ♦  /1X8HPFRI0D    9F12.5*/ 
1  lX8HrRFQ      9F12.5»/9H0AVG  FRFQ  lFip.e;  *SX5H0MEf!A  F12.5* 

1  3X9H0  XX  OV  C»E12.3) 

6240  IF(JK.GE.11)G0  TO  99Q99 

??9  00  2M  1  =  1. N 

7?T0  P(I»1 )aP(I.NJ  LESS  1 ) SU < I 1 1 ) =U ( 1 1 NJ  LESS  1) 

7231  P(I»?)=P(I»NJ)  $U(I«?)=U(I.NJ) 

231  CONTINUE 

232  GO  Tn  60 

92?n  F0RMaT(1X.I4»1HP  .16F8.0) 

R2?1  FORMAT(5X.iHU  ,i6F8.j) 

92??  FORMAT rSX»lHP  .16F8.C) 

999Q9  EmO   FrEOEnCY 
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